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We present a high-statistics calculation of nucleon electromagnetic form factors in Nf =2 + 1 
lattice QCD using domain wall quarks on fine lattices, to attain a new level of precision in systematic 
and statistical errors. Our calculations use 32 3 x 64 lattices with lattice spacing a = 0.084 fm for pion 
masses of 297, 355, and 403 MeV, and we perform an overdetermined analysis using on the order of 
3600 to 7000 measurements to calculate nucleon electric and magnetic form factors up to Q 2 ~ 1.05 
GeV 2 . Results are shown to be consistent with those obtained using valence domain wall quarks with 
improved staggered sea quarks, and using coarse domain wall lattices. We determine the isovector 
Dirac radius r", Pauli radius r^ and anomalous magnetic moment k„. We also determine connected 
contributions to the corresponding isoscalar observables. We extrapolate these observables to the 
physical pion mass using two different formulations of two-flavor chiral effective field theory at one 
loop: the heavy baryon Small Scale Expansion (SSE) and covariant baryon chiral perturbation 
theory. The isovector results and the connected contributions to the isoscalar results are compared 
with experiment, and the need for calculations at smaller pion masses is discussed. 
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I. INTRODUCTION 

Electromagnetic form factors characterize fundamental aspects of the structure of protons and neutrons, in particular 
they specify the spatial distribution of charge and magnetization. For non-relativistic systems the electric and magnetic 
form factors would just be Fourier transforms of the charge and current densities. At each Q 2 , the Sachs form factors 
Ge{Q 2 ) and Gm(Q 2 ) may be regarded as three dimensional Fourier transforms of charge and magnetization densities 
defined in the corresponding Brcit frame. A probabilistic interpretation of the Dirac and Pauli form factors Fi(Q 2 ) 
and F2(Q 2 ) can be obtained from a two dimensional Fourier transformation to impact parameter space in the infinite 
momentum frame [l|, HJ. At high momentum transfer, the elastic form factor specifies the amplitude for a single quark 
in the nucleon to absorb a very large momentum kick and share it with the other constituents in such a way that the 
nucleon remains in its ground state instead of being excited. It thus describes the onset of scaling and the scale at 
which quark counting rules become applicable, which is an unresolved theoretical question in nonperturbative QCD. 
The combination of precision experimental measurements and crisp theoretical interpretation renders elastic nucleon 
form factors particularly significant. Given the constantly improving experimental measurements of form factors and 
their fundamental significance, it is an important challenge for lattice QCD to calculate them accurately from first 
principles. 

The nucleon Dirac and Pauli form factors, Fi(Q 2 ) and F2(Q 2 ) respectively, are defined as follows for each quark 
flavor (/): 



P',S'\V{ l f) \P,S) = U(P',S') 



^F[ t \Q 2 )+ l ^^ F F^(Q 2 ) 



N 



U(P,S), Vf f) = ^ f)1 ^(f) , 
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where P, P' are the initial and final nucleon momenta, S 1 , S' are the corresponding spin vectors, the momentum 
transfer is q = P' — P with Q 2 = —q 2 > 0, and Mjy is the nucleon mass. The Sachs form factors Ge(Q 2 ) and Gm(Q 2 ) 
arc defined by: 

G E {Q 2 ) = Fl {Q 2 )-j0-^F 2 {Q 2 ) (2) 
G M (Q 2 ) = Fi(Q 2 ) + F 2 (Q 2 ) ■ (3) 

Finally, it is useful to define isoscalar and isovector form factors as the sum and difference of proton and neutron form 
factors as follows: 

FUQ 2 ) = Fl 2 {Q 2 ) - F" 2 (Q 2 ) = Fl 2 {Q 2 ) Fl d 2 (Q 2 ) = F^ 2 d (Q 2 ), (4) 
FUQ 2 ) = Ff. 2 {Q 2 ) + F{) 2 (Q 2 ) = I (F« 2 (Q 2 ) + F* 2 (Q 2 )) = \f$\Q% (5) 

where Ff' n are the form factors of the electromagnetic current in a proton and a neutron, respectively: 

Kro.p = \^fu - l*fd > Kro.n = + . (6) 

Although proton and neutron form factors contain both connected diagrams, calculated in this work, and disconnected 
diagrams, which are currently omitted, the disconnected diagrams do not contribute to the isovector form factors F?. 
Hence, we will devote particular attention in this work to the isovector form factors. 

Precise experimental measurements of the set of all four nucleon form factors remains challenging, and the field is 
marked both by significant recent developments and open questions. Although the most straightforward measurement 
is Fi(Q 2 ) for the proton, the slope at very small values of Q 2 remains controversial. Phenomenological fits to 
experimental form factors 0, 13 appear to be inconsistent with analyses based on dispersion theory 0,0.0. 

with 

phenomenological fits yielding larger Dirac radii. Hence, a new generation of precision measurements of form factors at 
low momentum transfer is currently being undertaken at Mainz Q . Spin polarization experiments 0, [ToL ITU IT2I HH 
yielded results for F 2 (Q 2 ) significantly different from traditional measurements based on Roscnbluth separation, 
and there is a consensus that two-photon exchange processes contribute much more strongly to the backward cross 
section used in Rosenbluth separation than to polarization transfer Q . However, there are not yet precise theoretical 
calculations of two photon exchange that fully resolve the discrepancy between the two experimental methods, and 
hence experiments usin g po sitron scattering, for which the relative contribution of the two-photon term changes sign, 
arc being prepared (Til Il5j|. Neutron form factors arc more uncertain than proton form factors because of the need 
to know the nuclear wave function to go from experimental scattering results from deuterium or 3 He to a statement 
about the neutron form factor. Over the years, nuclear models and theoretical calculations have been refined, but it 
is still a challenge to provide a definitive estimate of the uncertainty in the claimed neutron form factors extracted 
from nuclear targets. Given the level of precision to which we aspire in lattice calculations, systematic uncertainties 
in isovector and isoscalar form factors are not necessarily negligible. In the future when lattice calculations reliably 
include precise calculations of disconnected contributions, it may well be that lattice calculations play a role in guiding 
the resolution of some of these experimental questions. 

Electromagnetic form factors have now been calculated in lattice QCD using a variety of actions. Quenched 
calculations of form factors have used both Wilson [T(| [Tt} and domain wall [TU fermion actions, and additional 
quenched calculations have addressed magnetic moments and root-mean-squared (rms) radii Il9l. |20H . Dynamical 



calculations with two flavors have used Wilson |17( , clover improved Wilson |2l| , twisted mass |22j, |23j and domain 
wall |24| actions. Extensive 2+1 flavor calculations have been performed with a mixed action, which combines domain 
wall valence quarks and improved staggered sea quarks |25l. I26ll27j . using the same methodology as in the present work, 
and comparisons will be made to assess the consistency of the full domain wall and mixed action results. Dynamical 
domain wall results with 2+1 flavors on coarse lattices with a = 0.114 fm have recently been reported [2g, [H, |29| . 
and initial results from the present work on fine lattices with a = 0.084 fm were presented in Ref. [3C| . 

The goal of this work is to achieve a new level of precision in calculating form factors from first principles in lattice 
QCD. Hence, we have done everything feasible within the constraints of our computational resources to reduce both 
statistical and systematic errors. Since this involves a number of new developments, we describe our methodology, 
innovations, and tests in detail. Because the signal to noise for baryon observables degrades with increasing Euclidean 
time t as e -( M «- 3 / 2m f )* ; we have obtained high statistics using from 3688 to 7064 measurements of operators at a 
given mass by performing 8 measurements per lattice and have verified their statistical independence. The source-sink 
separation distance is a crucial issue, since an excessively large distance degrades the statistical accuracy whereas too 
small a distance introduces systematic errors from the contributions of excited states. We present a quantitative 
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analysis of the contributions of excited states, and using this analysis, provide compelling numerical evidence that 
with our choice, which has been questioned in the literature p9j . excited state contributions are negligible in our 
present work. Our overdetermined analysis of form factors provides a general framework for optimizing the precision 
of our lattice calculations by combining measurements of as many distinct nucleon matrix elements involving the form 
factors at the same Q 2 value as practical. We also describe how we choose which contributions to include, and treat 
error correlations. We compare domain wall calculations on fine lattices at three masses with a calculation on a coarse 
lattice at one mass, and present evidence that the 0(a 2 ) corrections are indeed small. We also compare our results 
with mixed action results, showing essential consistency between mixed and domain wall actions and emphasizing the 
small size of finite volume corrections to calculations on a 2.5 fm lattice at = 350 MeV that have been calculated 
to high precision with the mixed action. We perform chiral extrapolations of the Dirac and Pauli mean squared radii, 
(7y s ) 2 and (V^) 2 > respectively, and of the anomalous magnetic moments K„ jS . We use two different formulations of 
ST/ (2) chiral effective field theory: the heavy baryon Small Scale Expansion (SSE) which includes explicit A (1232) 



degrees of freedom 31 and covariant baryon chiral perturbation theory (CBChPT) without an explicit A (1232) in 
the Ji?-schcme 32|, 33 j, which represents a variant of infrared regularization [111 1 . We explore the degree to which 



the relevant low-energy constants can be determined in the range of masses we consider and the variation of the 
extrapolated results in both schemes. We conclude that with the new level of precision we achieve, it is necessary to 
extend the lattice calculations to substantially lower masses to make contact with the regime of applicability of chiral 
effective theory and possibly reach agreement with experiment. 

The remainder of this paper is organized as follows. In Sect.[ITl we present a detailed description of our methodology, 
including setting the scale, computation of nucleon matrix elements and coherent sink technique, optimization of 
sources, treatment of error correlations and constraints in the overdetermined analysis, and a check of the independence 
of multiple measurements per configuration. Sect. Hill presents the results of our lattice calculations for isovector form 
factors, including phcnomcnological fits to the momentum transfer dependence and determination of the Dirac radius 
(r") , the Pauli radius (r^) , and the anomalous magnetic moment k v . Comparisons are made with domain wall 
calculations on a coarse lattice and with mixed-action calculations using valence domain wall valence quarks and 
improved staggered sea quarks. We also present the chirally extrapolated values of (r\) 2 , (r%) 2 , and K v to the 
physical pion mass using the SSE and covariant chiral effective field theories and compare them with experiment. 
Corresponding results for isoscalar form factors are presented in Sect. HVl Systematic errors are discussed in Sect. [Vj 
results are compared with other work in Sect. IVI1 and conclusions and opportunities for further understanding of 
nucleon form factors are discussed in the final Sect. IViTl Selected numerical results are tabulated in Appendix lAl and 
the optimized sources are described in Appendix [B] 



II. LATTICE METHODOLOGY 
A. Dynamical domain wall ensembles 

In our calculations, we analyze gauge configurations generated by the RBC and UKQCD collaborations [I?) with the 
Iwasaki gauge action and Nf = 2 + 1 flavors of dynamical domain wall fcrmions. The gauge configuration ensembles 
are summarized in Tab. [I] We obtain the relevant physical results from three fine lattice ensembles with lattice spacing 
a = 0.084 fm. We use one coarse lattice ensemble with known lattice spacing a = 0.114fm [U to set the scale on the 
fine lattices and control the systematic errors due to discretization. 

In our analysis, we use only a unitary fermion action, where the sea and valence fermion actions and masses are 
exactly the same. The extent of the fifth dimension is chosen to be L s = 16, which keeps the residual mass m rcs 
smaller than the bare quark masses for all ensembles. 

In order to maximize the signal to noise ratio and suppress excited state contamination, we carefully optimize the 
quark propagator sources. We use Wuppertal smearing of quark sources combined with APE smearing of the source 
gauge fields to reach the maximum overlap of the lattice nucleon operators with the nucleon ground state and reduce 
its fluctuation. The details of optimization and the source parameters we use are given in Appendix iBl 

To increase statistics, we perform eight measurements of nucleon correlation functions on each gauge configuration. 
To do so, we compute four forward quark propagators and construct nucleon and antinucleon correlators advancing 
in the positive and negative time directions, respectively. The data for antinucleons are transformed according to 
the reflection symmetry and combined with the data for nucleons into a single data set. We save computing time by 



For recent work on chiral extrapolations of nucleoli magnetic form factors and octet-baryon charge radii in heavy baryon ChPT with 
finite range regularization, we refer the reader to Refs. 1 3 5|. 1 3 611 . 
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Table I: Gauge configuration ensembles used for our analysis, with one coarse and three fine lattice spacings. These configura- 
tions were generated by the RBC and UKQCD [33] collaborations. The coarse lattice spacing was determined in [37j], and the 
fine lattice spacing is determined in Sect. Ill Bl Measurement count includes a factor of 8 for each gauge configuration. Note 
that for m n , F n , m'^ the measurement count is the number of configurations multiplied by 4 instead of 8. 



Li x L t a [fm] T 


# 


ami/amh om! ves x 10 a am, [MeV] 


clF.* [MeV] 


aM N M N [MeV] 


24 3 x 64 0.114 9 


3208 


0.005/0.04 3.15(1) 


0.1901(3) 329(5) 


0.06100(11) 105.5(1.7) 


0.657(4) 1136(20) 


32 a x 64 0.084 12 
32 3 x 64 0.084 12 
32 3 x 64 0.084 12 


4928 
7064 
4224 


0.004/0.03 0.665(3) 
0.006/0.03 0.663(2) 
0.008/0.03 0.668(3) 


0.1268(3) 297(5) 
0.1519(3) 355(6) 
0.1724(3) 403(7) 


0.04400(15) 102.9(1.8) 
0.04571(09) 107.0(1.8) 
0.04755(18) 111.3(2.0) 


0.474(4) 1109(21) 
0.501(2) 1172(21) 
0.522(2) 1221(21) 



using the "coherent" backward propagator technique, in which we compute only a sum of four backward propagators 
for four separate sequential sources with the same hadron type, flavor and sink momentum. To check for possible 
systematic effects, we recalculate the nucleon three-point functions using independent backward propagators and 
larger source-sink separation on a subset of our lightest pion ensemble, and the extracted form factors (see Fig. |20|) 
show no significant deviation from the method we use. Since lattice data may be autocorrelated, we block all the 
measurements on the two consecutive gauge configurations, and also check that the measurements we get are indeed 
independent by increasing the block size to include eight consecutive configurations (see Fig. [2]). 



B. Pion mass, decay constant and setting the fine lattice scale 

So far, the scale has been set only for the coarse lattice ensembles [37} . In order to set the scale for the fine lattice 
ensembles, we compare the lattice values for the pion decay constant (aF„) on coarse and fine lattices at the same 
value of the dimcnsionless ratio (m 7r /i r 7r ) 2 ignoring possible finite lattice spacing effects in the pion decay constant 
F n . 

First, we compute the pion mass, the pion decay constant and the local axial current renormalization constant 2 
from fits to the pseudoscalar density and axial current correlators using the PCAC relation [39j : 

(A (t,p = 0)Jb(0)) = (e-«-« - e -M^-0) x 9 J^j. , x Z^Z^, (7) 
(J 5g (t,P = 0)J 5 (0)) = f e -m,t + e -m A L t - t) \ x F >* x m' tes Z-^ (8) 

\ / <±(j/lj -t Tn tes ) 

(J 5 (t,p= 0)J B (0)> = (e-^' + e-^^-'A x f * m * x Z£, (9) 
V / 4(m; + m[ cs y 

(J 5 (t,p= 0)Jb(0)) = (e-^' + e-"^-')) x - f * m * x Z" 2 , (10) 
\ / t(Trii -\- Tn rcs ) 

where Aq is the local axial charge, J$ q is the fifth dimension mid-point pseudoscalar density and J5 (J5) is the 
(smeared) pseudoscalar density. The pion decay constant convention is such that 

F ph ys = Q2 4 ± 3 McV ^ ^ 

We choose the range of t to be [12 : 52] to exclude any excited state contaminaions. We define the smearing 
renormalization constant Z sm from the plateau (Js(t) Js(0))/(J5(t) Js(0)) and the local axial current renormalization 
constant Za from the ratio of (Ao(t + 1/2)^(0)) and (An(i) Jr(0) ) ap propriately averaged to suppress 0(a) effects due 
to a 1 2 displacement of the conserved axial current Ao(t +1/2) [391 ]. The results for am^. aF^ and am' les are shown 
in Tab. [TJ The error bars reflect both the statistical error and the systematic error due to different fitting ranges. 
Second, we fit and F„ at three values of the light quark mass using 0{p A ) SU(2) chiral perturbation theory 



2 In this paper, we assume that the renormalization constant Za of the (partially) conserved domain wall axial current is equal to 
its L s — ► 00 value of one, and note that the finite L s deviation has been estimated in l37t |38H to give |Z_4 — 1| < 1%. The values of i*V 
that we compute are, in fact, F^/Za- 
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mm 

»X = a 2 X {l + 5 (0 + log } , (12) 

aF « = aF I 1 + (5^(0 - rn^Ff ^ ^ } ' ( 13 ) 

where et 2 % = 2(aB) ■ a(m; + mj. es ), ^^(a" 1 ) are the next-to-leading order (NLO) low-energy constants (LECs) at 
the scale A = a -1 , and the fit variables are (aF), (aB) and 1% 4 . However, the fit is not satisfactory in terms of x 2 : 
for two degrees of freedom, we get \ 2 w 7, with its probability to be this or higher being < 3%. This is the first 
of many indications that chiral perturbation theory, at the order we can use, is not accurate in the range of masses 
we are considering. Hence, the LEC's are not precisely determined although, as noted below, we obtain an adequate 
interpolation to set the scale. 
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Figure 1: One-loop SU(2) ChPT interpolation of the fine lattice values of F n and m w . The point with abscissa 9.71 = 
(m^/F^-)c 0arso was obtained by interpolating (aF n ) linearly in (m,/F,) 2 . 

The NLO LECs /g 4 from our fit can be converted to the scale-independent parameters l^^ [4(|. At the physical 
pion mass we obtain 

h = 3.08(11), h = 4.24(4). (14) 

Our result for Is i s in agreement with the crude estimate I3 = 2.9 ± 2.4 [4(j an d with the lattice determination 
[3 = 3.0(5)(1) m usm g Nf = 2 dynamical Wilson fermions but disagrees with l 3 = 3.42(8)(10) (the errors are 
statistical and systematic due to residual lattice artifacts) from the ETM collaboration [43[ . This discrepancy could 
could arise from the difference between Nf = 2 and Nf = 3 flavors of dynamical fermions. Furthermore, chiral 
symmetry implies that I4 determines the slope of the scalar form factor of the pion. In their seminal paper, Gasser 
and Leutwyler obtain I4 = 4.3 ± 0.9 [4(j. This estimate has been sharpened in (44|: I4 = 4.4 ± 0.2, which agrees with 
the value l 4 = 4.4 ± 0.3 obtained by Bijnens et al. [H- The E ™ collaboration result [H[ is l 4 = 4.59(4) (2). 

The resulting interpolated functional dependence of (aF„) on (m,/^) 2 is shown in Fig.[TJ For simplicity, we also 
estimated (ai*^)!* at {m^ / F^) 2 \ arsc and its error by linear interpolation in (m^jF^) 2 between the two lightest pion 
masses. The comparison in Fig. |¥jshows no difference between these two approaches. We also interpolated the lattice 
value of the nucleon mass (aNpj), and obtained the ratios 

{aF ^* - 0.735(2), = 0.742(5). (15) 



coarse v " (aM N )\ c 

Although these ratios are barely consistent within errors, their discrepancy is irrelevant to the fine scale determination 
as long as the fractional error in the coarse lattice scale a coa rse = 0.1141(18) fm dominates. We obtain the value for 
the fine lattice scale 

Ofine = 0.0840(14) fm, a^ c = 2.34(4) GeV. (16) 
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C. Extraction of nucleon matrix elements 

In order to calculate nucleon matrix elements, we compute the three-point polarized nucleon correlators involving 
the vector current, along with the two-point correlators [251 ] : 

C 2pt (t,P) = ^e- ii3 -^(r pol ) a/3 (^(f,t)iV a (0,0)), (17) 

x a/3 

CZ (r, T; P, P>) = J2 e -i?'-*+<?'-PyitJ2 (r pol ) Q/3 (Np(x , T)W% r)N a (0, 0)) (18) 



where Np,N a are the lattice nucleon operators, (f2 |iV a (a;)| P, a) = ^/Z{P)UZ '{P)e' tPx , with Z{P) parameterizing 
the overlap with the nucleon ground state, (r po i) Q ^ = 1 ~t ; 74 1 ~'^ 375 is the spin and parity projection matrix 3 , and 
yan _ q^Hjf'q j s ^j^g vector current operator, where t a denotes an isospin generator. In the transfer matrix formalism, 
these correlators take the form 

Z(P)e~ Et 

C 2p t(t, P) = Tr [r po i {if + M N )] + excited states, (19) 

C%£(t, T; P, P') = ^ Z{P,) ' Z ^) e 2 ^ iT T) [ r poi {if' + Mn) T» (P' : P) {if + Mjv)] + excited states, (20) 

where E and E' are the ground state energies of the initial and hnal nucleon states and T M (P', P) is the electromagnetic 
vertex function defined below in Eq. ([24]) . Excited state contributions have generally similar forms with different Z- 
factors, vertices and higher energies E cxc > E. The systematic effects related to them will be discussed in Sect. IV Al 



Table II: Momentum combinations used to extract the form factors (only one representative of in/out momenta is given). 
Approximate Q 2 values are given for the lightest Mn = 1109 MeV. 



# 


(out | in) 


Q z [GeV*] 


1 


(0,0,010,0,0), (-1,0,01-1,0,0) 


0.0 


2 


(0,0,0|1,0,0), (-1,0,0|0,0,0) 


0.203 


3 


(-1,0,01-1,0,1) 


0.204 


4 


(0,0,0|1,1,0) 


0.391 


5 


(-1,0,01-1,1,1) 


0.395 


6 


(-1,0,0|0,0,1) 


0.422 


7 


(0,0,0|1,1,1) 


0.568 


8 


(-1,0,0|0,1,1) 


0.626 


9 


(-1,0,0|1,0,0) 


0.844 


10 


(-1,0,0|1,1,0) 


1.048 



In order to extract the combinations of matrix elements (P 1 , S' \V^\ P, S) = U{P', S')T»{P', P)U{P, S), we combine 
the lattice nucleon correlators (fl~9ll20|) into the usual ratio of 3- and 2-point correlation functions, which we find useful 
to write in a convenient and illuminating new form as follows. First, we define two ratios, a normalization ratio, i?j\r, 
and an asymmetry ratio, Ra, 

Rn = , 3p , (21) 

y/C2 pt (T,P)C2 P t(T,P')' 

j c 2pt {T-T 7 P)C 2pt {r 7 P') 

~ \jc 2pt (T-T,p')c 2pt (r,py [2Z) 

The physical matrix element is then given by the product: 



3 In this subsection, we use Euclidean 7-matriccs, (7 M )t = 7 M , {7 M ,7^} = 2<5 M ". 
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R v " ee R N R A = C ^ (T ' T; P ' P ' ] I ° 2pt {T ~ T > P) ° 2pt (T ' Pl) 



y/C 2pt (T,P)C 2pt (T,P')y C2 P t(T - r, P')C 2pt (r, P) 
,00 Es,s< {U(P,S)T pol U(P',S')) ■ (P',S'\V,\P,S) 



^/2E{E + M N ) ■ 2E'{E' + M N ) 

The normalization ratio, Rn, has the property that all the lattice-dependent overlap factors Z for the ground 
state cancel out, which motivates its name, and yields the full result in the case of forward matrix elements P = P' . 
The asymmetry ratio, Ra, compensates the asymmetric exponential t dependence of the three-point correlator, 
which motivates its name. In the absence of excited states, it would be equal to exp [— [E 1 — E)(t — T/2)\ and in 
the forward case, P' = P, this ratio is trivial and equal to one. In the general case, P' 7^ P, this ratio is still 
identically one in the center of the plateau, r = T/2, and possesses the following symmetry around the plateau center: 
Ra{T~t) = 1/R a (t). 

The limit T — > 00 should be taken to get rid of the excited state contamination. In practice, this requires adopting 
a value of source-sink separation T large enough so that the excited state contributions to Eq. (|23|) arc negligible 
compared to the other sources of errors. We will explicitly explore the contributions of excited states to R v> *m 
Section [Vj where the decomposition into the product RnRa will prove extremely useful. 

In order to obtain the most precise information on the form factors, we constrain the in- and out- lattice nucleon 
momenta to have components 0,±1. Higher momentum components are subject to stronger finite lattice spacing 
effects, i.e., discretization errors and dispersion relation deviations from the continuum expression. There is also an 
indication (see Sect. IV K\ that such states have larger excited state contaminations. 



D. Overdetermined analysis of form factors 

In Minkowski space, the nucleon electromagnetic vertex T^(P',P) in Eq. (|2"0)) is parameterized with two form 
factors: 

T»(P\P)=F 1 (Q 2 )^+F 2 (Q 2 ) 1 -^^, q = P'-P, Q 2 = -q 2 . (24) 

Transforming the above expression to Euclidean space and substituting it into Eq. (f2T)|) and then Eq. (|2"5)) and 
neglecting the excited states, we obtain an overdetermined system of equations for the form factors F lt2 {Q 2 ) at each 
fixed value of Q 2 : 

A ai Fi(Q 2 ) = R^, a = 1,2,... (25) 

where we use a summation convention over i = 1, 2 and a is a composite index specifying the current component and 
the initial and final momenta of a given matrix element (for fixed Q 2 ), which will be discussed below. The r.h.s. of 
Eq. |25|) is evaluated using Eq. |23|) with computed lattice correlators. 

We find the solution of the overdetermined system from a linear fit, which minimizes the functional 

T = J2 (AaiFi - R a ) Cj {A P] F 3 - R P ) , (26) 

a/3 

where C a p is the covariancc matrix of R a averages, C a p = jf^j (((R a Rp)) - ((Ra})((Rp)}), with the double brackets 
denoting an ensemble average. Using the covariance matrix is crucial as long as the correlation functions prove to be 
correlated. 

Since the covariance matrix may be ill-determined, it can introduce uncontrollable errors into the extracted form 
factors. In general, a covariance matrix is notoriously difficult to reliably estimate in a statistical analysis. To make 
sure the linear fitting gives a correct result, we repeat the analysis with only the diagonal elements of the covariance 
matrix C QQ , which is equivalent to an uncorrelated linear fit. The comparison of these two schemes is presented in 
Fig. [2 We find that the form factors from an uncorrelated fit are consistent with the correlated fit results. 

The overdetermined system (|25[) contains a subclass of equations which have an exactly zero l.h.s.: A a i = 0, i — 1, 2. 
The measured lattice value of a right-hand side R a is not required to be zero, and may be correlated with other matrix 
elements. In an uncorrelated fit, such equations decouple and do not contribute to the solution. In contrast, the 
outcome of a correlated fit depends on such values, thus potentially better utilizing the input from lattice calculations. 
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In addition, by fitting the equations with a vanishing l.h.s., we check the symmetries of the electromagnetic vertex 
statistically. Fig. [5] also shows the agreement of the full overdetermined system solution and the system without 
zero l.h.s. equations, confirming the consistency of our analysis. 

The dimension of the overdetermined system may grow large, especially when many momentum combinations are 
included. For example, the most precise point for Q 2 > corresponds to the matrix element (0,0,0 |V^ M (0)| 1,0,0). 
All components, together with spatial rotations and reflections give 48 equations, only 16 of which are non-zero. It 
is useful to combine all the nuclcon matrix elements for fixed Q 2 into equivalence classes based on spatial (rotational 
and reflection) symmetry. We adopted the following heuristic equivalence criteria 4 for three-point functions: 

• The momenta of the in- and out-states must be equivalent under the spatial symmetry. 



• The corresponding coefficients A a i in Eq. 



must be equal up to an overall sign. 



• The component of the current operator must be temporal or spatial and real or imaginary for both matrix 
elements being compared. 

Blocking the measurements in each equivalence class is advantageous for two reasons. First, this reduces the 
dimension of the system of equations (|25|) and the covariance matrix we need to estimate, and we note that blocking 
strongly correlated values improves the covariance matrix condition number. Second, as long as for the equivalent 
three-point functions we need spatially equivalent two-point functions to build the ratio in Eq. (|23|) . we can block the 
two-point functions separately before computing the ratio. This improves the method in Eq. (|23|) by reducing the 
fluctuations of the two-point functions in the denominator. 

To extract the final set of the form factors, we perform a correlated fit to the reduced (i.e., the system with no 
equations whose l.h.s. is zero) overdetermined system with blocked equivalent equations. 
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Figure 2: Comparison of the nucleon form factors extracted from the full overdetermined system, only non-zero equations, 
uncorrelated fit and averaged equivalence classes for — 297 MeV. Increased binning of data (eight successive configurations 
instead of two) shows no increase in estimation of statistical errors. Each form factor value is divided by the central value of 
the dipole fit. Tab. HT1 lists the momentum combinations corresponding to each index on the horizontal axis. 



III. ISOVECTOR FORM FACTORS 



In experiments, the proton and neutron electromagnetic form factors (see Eq. ([6])) are measured separately, and the 
isovector form factors (0| can be calculated by taking their difference. In lattice calculations, the Wick contractions of 
the quark fields in Eq. ([6]) with nucleon operators indicate that disconnected quark loops in the current insertion would 
be needed to calculate the proton and neutron form factors separately The calculation of the disconnected quark 



4 We have not classified the matrix elements according to the hypercubic lattice symmetry but instead use relations derived in the 
continuum. Thus these criteria may be thought of as numerical means to improve the condition number of the linear system we need 
to solve. 
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loops is numerically demanding and has not been included in current calculations. However, the disconnected loop 
contributions cancel (in the isospin limit) in the contraction of the difference of the proton and neutron electromagnetic 
currents in Eq. ([6|), which gives the matrix elements needed for the isovector form factors. We focus our discussion 
on the isovector form factors in this section. 

After presenting our lattice results for the isovector Dirac and Pauli form factors and the rms radii, we will compare 
chiral extrapolations using the SSE formulation and covariant baryon chiral perturbation theory 5 . Corresponding 
results for the connected contributions to the isoscalar form factors will be presented in Sect. IIV1 

A. Vector current renormalization 

The isovector Dirac form factor at zero momentum transfer, Fi~ d (0), gives the difference of the electric charges for 
the proton and neutron, which is 1. Since we can measure _F"~ (0) very accurately on the lattice, we use it to obtain 
the vector current renormalization constant, Zy, by setting 

Zv-Fi _d (0) = 1. (27) 

Since domain wall fermions have good chiral symmetry, in the chiral limit the vector current renormalization, Zy, and 
the axial vector current renormalization, Za, are expected to be the same up to 0(a 2 ) corrections. Za is measured 
by taking the ratio of the point-split five-dimensional conserved axial current to the local four-dimensional current 
(see Sect. Ill Bl and [IH). We show the results of Zy and Za in Tab. IIII1 Naive linear extrapolations in m\ to the 
chiral limit show that Zy and Za are consistent within errors, as is clearly shown in Fig. [3] 

Table III: Vector and axialvector current renormalization constants. The chiral limit values are obtained by linear extrapolations 
to mjj. = 0. 

[MeV] Zy Z A ~ 

297 0.7468(39) 0.745025(24) 
355 0.7479(22) 0.745207(18) 
403 0.7513(17) 0.745317(20) 
chiral limit 0.7397(74) 0.744700(55) 



0.755 1— r 




[GeV 2 ] 

Figure 3: Comparison of the vector and axial vector current renormalization constants. In the chiral limit, these two renor- 
malization constants agree within errors. The errors on all the Za points given in Tab. IIIII are too small to appear on the 
figure. 



5 For an analysis of nucleon electromagnetic form factors in baryon ChPT with standard infrared regularization, we refer to [4(j | . 
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In the following analysis, we renormalize the form factors by Zy as measured on the corresponding ensemble. That 
is, we use a mass-dependent renormalization condition. The mass dependence of the renormalization constants is very 
mild and consistent with the theoretically expected form Zy(g )(l + b v am q ) [47j |. 



B. Q 2 dependence 

As will be discussed in the following section, ChPT describes the Q 2 -dependence of the form factors for values of 
Q 2 much less than the chiral symmetry breaking scale (typically of the order of the nucleon mass). Lacking a model- 
independent functional form applicable in the largc-Q 2 region, we study the Q 2 dependence using the phenomenological 
dipolc or tripole formula. The Dirac form factor is fixed to 1 at Q 2 = under our renormalization scheme, and we 
use the following one-parameter dipole or tripole formula to describe the Q 2 dependence: 



Fi(Q 2 ) = ^2 (one-parameter dipole), (28) 

Fi(Q 2 ) = q2 (one-parameter tripole). (29) 

The Pauli form factor at Q 2 = 0, -F^O), cannot be measured on the lattice directly. We thus fit the data using the 
two-parameter dipole or tripole formula, 

-^MQ 2 ) = -772 (two-parameter dipole), (30) 

F2(Q 2 ) = — (two-parameter tripole). (31) 

We are interested in mean squared Dirac and Pauli radii, which are defined by the slope of the form factors at small Q 2 : 



Fi, 2 (Q 2 ) = Fi, a (0) 
and are related to the pole masses by 



l- l -{r l .2) 2 Q 2 + 0{Q i ) 
6 



(32) 



<r)' - (33) 



for the dipole fits, and 



(rf = —^-2, (34) 



for the tripole fits. 

Note that results at different Q 2 from the same ensemble may be highly correlated [26|, therefore we perform 
correlated least-x 2 fits to the data. We investigate the extent to which the dipole and tripole Ansiitze describe our 
data and the stability of the fits by varying the maximum Q 2 values included in the fits. 

In Tab. IXVTl we show the fit results for F™~ d {Q 2 ) using the one-parameter dipole and tripole formulae in Eqs. I]28p 
and (|29|) . Comparing the x 2 /dof for the dipole and tripole fits, we see that the dipole fits are slightly preferred when 
larger Q 2 values are included in the fits. However, the Dirac radii determined from both the dipole and tripole fits 
agree within errors. In general, the dipole form describes the data reasonably well throughout the whole Q 2 range for 
all but one ensemble, the = 355 McV ensemble, where, when Q 2 cutoff is larger than 0.3 GeV 2 , x 2 /dof becomes 
very large. This may be due to the fact that this ensemble has the most statistics, and we start to see the deviation 
from the phenomenological dipole formula. For the other two ensembles, we can see the general trend that when 
large Q 2 points are included in the fits, the x 2 /dof becomes slightly worse, while the fit parameters do not depend 
significantly on the choice of the Q 2 cutoff, indicating that the dipole fits are stable. 

We do the same comparison for F^ _d ((5 2 ) as shown in Tab. IXVIIl Judging from the % 2 /dof values, we do not see 
significant differences between the dipole and tripole fits. Since the Pauli form factor is not constrained at Q 2 = 0, 
including larger Q 2 in the fits does not seem to affect the quality of the fits significantly. The fit parameters ^(0) 
and Mo,t prove not to be affected as well. 



11 



As an example, we show the dipole fit curves with a Q 2 cutoff at 0.5, 0.7 and 1.1 GeV 2 for the m n = 297 McV 
ensemble in the top panel of Fig. 2J To show the quality of the fits more clearly, we plot the ratios of the form factor 
data to the dipole fit with the Q 2 cutoff at 0.5 GeV 2 in the bottom three panels of Fig. |U The error bands reflect 
the jackknife errors in the dipole fit parameters. We see that although the data included in the fits can be described 
reasonably well by the dipole formula with discrepancies that are generally within two to three standard deviations, 
the clear systematic tendency indicates that the dipole Ansatz is not a good description of the data over the whole 
momentum transfer region. In particular, for F{~ d , the precisely measured points in the region of 0.2 GeV 2 are 
systematically lower than the dipole fit, whereas at high Q 2 , the lattice data are systematically higher. For 
the high Q 2 lattice data are systematically lower than the dipole fit. This is consistent with the empirical fits to 
the experimental data in Refs. [H, 3, where the phenomenological corrections to the dipole form are negative in the 
region of 0.2 GeV 2 and positive at about 0.4 GeV 2 . For comparison, we also plot the dipole fits with Q 2 cutoffs at 
0.7 GeV 2 (dashed line) and 1.1 GeV 2 (dotted line) relative to the 0.5 GeV 2 dipole fit (solid line). The differences 
between different Q 2 cutoffs are small, indicating that the fits are stable. 

It is worth noting that the Dirac and Pauli radii, r\ and r%, and the anomalous magnetic moment, k v , are defined 
in the Q 2 = limit. We thus restrict the fits to the smallest Q 2 points possible to extract these quantities while still 
including enough data points to constrain the fits. For uniformity we choose to determine these quantities from the 
one-parameter dipole fits for F^~ d , and the two-parameter dipole fits for i 7 2 " _d , with a Q 2 cutoff at 0.5 GeV 2 . 

We also perform dipole fits to Ge{Q 2 ) and Gm(Q 2 ) to see how well the dipole Ansatz describes the data. We find 
that the dipole fits to G^T d and G\J d are qualitatively similar to F™~ and F^~ d ■ However, it appears that the fits 
are even more stable over the whole range of Q 2 than Dirac and Pauli form factors. This is indicated by little change 
in the ratio plots in Fig. [5] with different Q 2 cutoffs. 

Figure [5] shows a comparison of the lattice results for Gb at three different pion masses from the fine ensembles 
and one pion mass from the coarse ensemble with a phenomenological fit to the experimental data using the pa- 
rameterization in Ref. (48j (with no indication of the experimental errors). The solid curves are dipole fits to the 
form factor results with the Q 2 cutoff at 0.5 GeV 2 . As the pion mass decreases, the slope of the form factors at the 
small momentum transfer monotonically increases. The results from the coarse ensemble at m„ = 330 MeV is nicely 
surrounded by the results from the fine ensembles at = 297 and m w = 355 MeV, indicating that the effect of the 
finite lattice spacing error should be small. 



C. Chiral extrapolations 

1. Chiral extrapolations using C(e 3 ) small scale expansion 



To compare the lattice results for the nucleon form factors at finite momentum transfer with the experimental 
results, we need to do extrapolations for both the m„ and Q 2 dependence using baryon chiral perturbation theory. 
This combined dependence has been worked out both in SSE at leading one loop accuracy and in CBChPT up to 
NNLO order. In particular, the 0(e 3 ) expression for the isovector Dirac form factor F^~ d (Q 2 , m n ) has been derived 
in Ref. [4l| and is given by 



Fr d {Q\m^ = 1 + 



1 



(471^)2 



dx 



dx 



16 



A 2 c\ 



81 

2 



40 2 5 2 

Q [ 7^ C A ~ 3SU 



35a 



+ Q z x(l - x) 5ffi + 1 - 



32 

dx —c\A ^A 2 - ml logi?(77v) - \/ A 2 - m 2 log R(m) 



1 

3 
40 



3] log 





log 


r ~ 9 1 

Tfl 


) 








m l. 



(35) 



where 



A 



R(m) = — + d— -1, 



A 2 



2 , 



mi + Q x(l - x). 



(36) 
(37) 
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Figure 4: The top panel shows the lattice results for F^ 2 d (Q 2 ) at = 297 MeV along with the dipole fits with three different 



Q cutoffs. The bottom left three panels show the ratios of the lattice results for _F" to the dipole fits using Eq. (|28[) . and 



the bottom right three panels show the ratios of the lattice results for to the dipole fits using Eq. (|30[) . Only the solid 

data points are included in the fits with cutoff 0.5 GeV 2 , and the grey bands show the errors for these fits. The dashed and 
dotted lines show the ratios of dipole fits at cutoffs 0.7 GeV 2 and 1.1 GeV 2 relative to the fit at 0.5 GeV 2 . 
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Figure 5: The top panel shows the lattice results for G U E M (Q ) at = 297 MeV along with the dipole fits with three different 
Q 2 cutoffs. The bottom left three panels show the ratios of the lattice results for G^T d to the dipole fits using Eq. (I28|l . and 
the bottom right three panels show the ratios of the lattice results for G^~ d to the dipole fits using Eq. (|30p . Only the solid 
data points are included in the fits with cutoff 0.5 GeV 2 , and the grey bands show the errors for these fits. The dashed and 
dotted lines show the ratios of dipole fits at cutoffs 0.7 GeV 2 and 1.1 GeV 2 relative to the fit at 0.5 GeV 2 . 
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Figure 6: Lattice results for G^T at three pion masses from the fine ensembles and one pion mass from the coarse ensemble, 
compared with a phenomenological fit to the experimental data as parameterized in Ref. [4^]. The solid curves are the dipole 
fits to the form factor results with a cutoff at Q 2 — 0.5 GeV 2 . 



In the above expressions, F„ denotes the pion decay constant in the SU (2) chiral limit with the convention in Eq. (JTTJ) . 
Here qa is the nucleon axial charge in the SU(2) chiral limit, ca is the leading-order pion-nucleon-A coupling 6 , and A 
denotes the A (1232)-nucleon mass splitting in the SU{2) chiral limit. For more details on the effective Lagrangians 
and the definitions of the low-energy constants, we refer the reader to [4lj . 

To the same order, the expression for the isovector Pauli form factor, is also derived in [ijj and is given as 



where, to C(e 3 ), 
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(38) 



(39) 



In order to capture the most prominent (Dim 2 ,) corrections, Hemmert and Weise |49l ] proposed a modification of 
the standard SSE power counting to promote the leading term of the magnetic N — > A transition into the first order 
iVA effective Lagrangian. This leads to the following expression for /^(nv): 
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(40) 



The coupling c A corresponds to g^NA ln the notation of Ref. [3lll . 
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where cy is the leading magnetic photon- nucleon- A coupling in the chiral limit and k° denotes the anomalous magnetic 
moment in the SU(2) chiral limit. We will use this expression in our analysis. 

Our results for the form factor F2 are given in terms of a quark mass dependent "magneton" (see Eq. (f24|) ). which 
is not accounted for in SSE at the order at which we are working. Therefore, in order to compare Eq. (|40[) with our 
lattice data, we follow Refs. [H, [lTj and define n nmm measured relative to the physical magneton: 



IV1 N Mat _ N 



-K 



■F 2 (0). (41) 



\ flat \ flat 

m N m N 

We then identify Mn in the SSE expressions as the physical nucleon mass. In the following comparisons of our results 
with chiral perturbation theories, the normalized magnetic moment K" orm will be used throughout, and we drop the 
superscript "norm" unless there is an ambiguity. 

ChPT describes the Q 2 -dependence of the form factors for values of Q 2 much less than the chiral symmetry breaking 
scale (typically of the order of the nucleon mass) and Q 2 counts as a small quantity, of the order of m 2 . In fact, we 
have attempted simultaneous fits to both the and Q 2 dependences of F^~ d using the SSE formula in Eq. (|35[) . 
and found that the fits fail to describe data even with Q 2 < 0.4 GeV 2 (x 2 /dof ~ 10). This is consistent with the 
findings of Ref. [4l[, where the applicability of the C(e 3 ) SSE results for the isovector nucleon form factors at physical 
pion mass was found to be limited to Q 2 < 0.2 GeV 2 . Lacking a model-independent functional form applicable in 
the large-Q 2 region, we resort to studying the pion mass dependence of the mean squared Dirac radius, (r^) 2 , Pauli 
radius, (r^) 2 , and the anomalous magnetic moment, k v , as obtained from the dipole fits discussed in Sect. [HTBl We 
tabulate these values in Tab. ITVl 



Table IV: Results for the isovector Dirac and Pauli radii and anomalous magnetic moment from dipole fits with Q 2 < 0.5 GeV 2 . 



m n [MeV] (rlf [GeV~ 2 ] 


{rlf [GeV- 2 ] 


< orm • (rS) a [GeV~ 2 ] 


norm 
K v 


297 7.83(21) 
355 7.23(14) 
403 6.98(13) 


9.82(84) 
9.55(46) 
9.74(41) 


24.1(3.0) 
24.1(1.7) 
24.5(1.5) 


2.447(99) 
2.518(57) 
2.508(51) 


330 7.46(22) 


11.44(67) 


31.6(2.7) 


2.758(84) 



The 0(e 3 ) SSE formulae for (r\) 2 and (r%) 2 can be derived from Eqs. (|35|) and |38|) . respectively, and are given by 




(42) 



%k 2 F 2 sJW^ 



(43) 



Together with the expression for the anomalous magnetic moment in Eq. (|40|) , these formulae involve six low- 
energy constants: F n , A, ca, $a, and cy, as well as two counter terms: -B 10 (A) and El(X). Ideally we would like to 
determine all these constants from simultaneous fits to lattice results. However, this is not feasible with the limited 
number of measured observables and pion masses in the present calculation, and we thus fix some of the low-energy 
constants using their phenomenological values. We describe our choices for these values below. 

In Ref. Colangelo and Diirr analyze numerically the NNLO expression for the pion mass dependence of F^ 
[45I ] . They use available information from phenomenology to fix all low-energy constants but the chiral limit value of 
Fn, use the physical value (fTT|) and obtain 



F~ 



(86.2 ± 0.5) MeV. 



(44) 



77 I chiral limit 

In the absence of reliable chiral extrapolations of both nucleon and A (1232) masses (see the discussion in Ref. [5ll]) 7 . 
we identify the A-nucleon mass splitting in the chiral limit with its value at the physical m v . The position of the 



7 For an analysis of the quark mass dependence of nucleon and delta masses in the covariant SSE at order e 4 we refer to loll . 
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A (1232) resonance pole in the total center-of-mass energy plane has been determined from magnetic dipole and 
electric quadrupole amplitudes of pion photoproduction. According to the Particle Data Group average [53|, the 
A-polc position leads to M A = (1210 ± 1) McV and T A = (100 ± 2) MeV. If one instead defines the A (1232) 
mass and width by looking at the 90° irN phase shift in the spin-3/2 isospin-3/2 channel, the PDG averages give 
M A = (1232 ± 1) MeV and T A = (118 ± 2) MeV. With M N = 939 MeV, one obtains, respectively, 

A = (271 ± 1) McV, (45) 



A = (293 ± 1) McV. (46) 

The A (1232) decays strongly to a nucleon and a pion with almost 100% branching fraction. From the PDG values 
of masses and widths [53| and from 

= l2lT p A 2 MA (K - "4) 3/2 (Ma + M N - E„), (47) 



where 



one obtains, respectively, 



2 n r2 , m 2 



Mi - M 



2M A 



\c A \ = 1.50... 1.55 if T = (100 ± 2) MeV and A = (271 ± 1) MeV; (49) 
\c A \ = 1.43... 1.47 if T = (118 ±2) MeV and A = (293 ± 1) MeV. (50) 

Calculating the strong decay width of A (1232) to leading order in (non-relativistic) SSE kinematics, one obtains 

r A ^ = ^(A 2 -m^/ 2 . (51) 

We note that this expression corresponds to the leading term in a 1/Mn expansion of the result given in Eq. (|47|). 
which utilizes the full covariant kinematics. Using the ranges of masses and decay widths mentioned above, this 
expression yields the lower values 

\c A \ = 1.11... 1.14 if T = (100 ± 2) MeV and A = (271 ± 1) MeV; (52) 
\c A \ = 1.04... 1.07 if T = (118 ±2) McV and A = (293 ± 1) McV. (53) 

Furthermore, SU (4) spin-flavor quark symmetry gives c A = 3 qa/ 0V2) = 1.34. 

Chiral extrapolations of different sets of lattice results [HI, HH, [56|, [57| based on SSE at leading-one-loop accuracy 
lead to a chiral limit value for g A of about 1.2. From the relativistic tree-level analysis of the process of pion 
photoproduction at threshold 7p — * n°p, one obtains [H,[5^] (for g^NA = 1-5) 

c v = (-2.5 ± 0.4) GeV" 1 . (54) 

As specified above, at the order C(e 3 ), all the couplings in Eqs. (|3"STl4"3"|) arc meant to be taken in the chiral limit. 
Replacing them with the corresponding quantities at the physical point amounts to the inclusion of higher-order 
effects. As long as the deviation between the values in the chiral limit and at the physical point is small, one expects 
such a replacement to yield little effect. To test this statement, in some cases we have performed the chiral fits using 
both the physical values and the chiral limit values for the low-energy constants and found no significant differences. 
In the following we will only present results obtained using the chiral limit values as inputs, which are summarized in 
Tab.lVl 



Table V: Input values for the low-energy constants in the fits. 

g A [GeV] A [GeV] 
1.2 0.0862 0.293 
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Among the low-energy constants discussed above, ca and cy are the two least known. In addition, we have little 
knowledge of the counter-terms, -B[ (A) and -E[(A), as well as the anomalous magnetic moment in the chiral limit, 
Ky, from phenomenology. Lattice calculations in the chiral regime have the potential to constrain these parameters 
to unprecedented accuracy. Our attempt here is to check the consistency of our data with the predictions of chiral 
effective field theories, to estimate the range of applicability of the ChPT formulas, and to determine these low-energy 
constants when the formulas are applicable. Since ca appears in the formulas for {r\) 2 , (r^) 2 and k v , a simultaneous 
fit to all these three quantities would give a better constraint for the value of ca- However, we have only three data 
points for each of these quantities, and k v alone has four parameters, three of which (cy, £l(A) and reJJ) are n °t 
constrained by any other quantity. Thus the quark-mass dependence of k v cannot be used to constrain ca- Therefore 
we choose to fit simultaneously 8 only (r") 2 and k v ■ (r%) 2 to determine ca and S[ (A), and then use the resulting ca 
as an input for the fit to n v . This way the three free parameters in n v are exactly specified by the three data points. 

We present the resulting x 2 /dof and fit parameters normalized at scale A = 600 MeV in the first row of Tab.[VT]and 
plot the fit curves as the solid lines in Fig. [71 As indicated by a x 2 /dof of 17, the simultaneous fit to (r\) 2 and k v -(r%) 2 
does not describe the data. The problem is that our results for (r\) 2 and K v ■ {r^) 2 favor different values for ca- In 
fact, an independent fit to (r") 2 yields ca = 1.98(7), while an independent fit to K v ■ (r%) 2 gives ca = 1.39(10). The 
tension between these two quantities results in the large x 2 /dof in the simultaneous fit, indicating that the formulae 
given in Eqs. (|42)) and (|43)) do not describe our data consistently. As we can see from Fig. |7(b)| the solid fit curve 
lies systematically higher than the data points, which then motivates us to add the (9(m°) correction to the leading 
one-loop result of Eq. (03]) (the so-called "core" contribution in Ref. [l6|) to k v ■ (r%) 2 , such that 



K v (m v ) ■ (r : 



v\2 



JY 



9aMn 



log 




2AM N C. 



(55) 



With this modification, the simultaneous fit to {r\) 2 and k v ■ (r%) 2 , now using Eqs. (|42)l and (|55"1) . appears to describe 
the average value of the data much better, but still not the pion mass dependence. We show the results in the second 
row of Tab. IVI| and the fit curves (dashed lines) in Fig. [71 The fit describes {r\) 2 very well, but ca turns out to be 
larger than the range discussed earlier, which, not surprisingly, gives rise to a smaller extrapolated value for {r\) 2 
than the exp eriments. Our new DWF data extend the trend of the weak pion mass dependence in (r%) 2 observed 
in Refs. [l6l . H3] now down into the range of pion masses ~ 300 MeV. The appearance of such a "plateau-like" 
behavior down to such light pion masses, which was also observed in Ref. [2||, is surprising. The leading one- loop SSE 
formulae (]43U55[) for this radius cannot accomodate such a behavior, with or without the inclusion of the higher-order 
"core" term. 

Using ca determined from the above fits either with or without the constant term in Eq. (|4"3"]) to {r\) 2 and k v - (r^) 2 , 
we fit k v to Eq. (|40[) with three unknown parameters, cy and E\{X). The results are shown in Tab. lVll The value 
for cy from our fit turns out to have a different sign from that determined in (58l . mentioned earlier. This is not 
surprising given that we only have three data points, which have little or no pion mass dependence. We do not have 
the freedom to check the consistency of the fit, and we do not expect to obtain a reliable estimation for cy, which, 
judging from Eq. ([4"0)) , is very sensitive to the curvature of the data. 



Table VI: Fit parameters from the fits to the isovector Dirac radius (»"i) 2 , Pauli radius (rj) 2 and the anomalous magnetic 
moment k v - Details of the fit procedures are described in the text. We have set the scale to A = 600 MeV. 





X Vdof 


CA 


c v [GeV" J 


K v 


Bio (A) 


E{(\) [GeV~ a ] 


C [GeV -3 ] 


no constant term 


17.0(4.0) 


1.54(6) 


8.7(5.8) 


4.13(95) 


1.20(17) 


-4.67(42) 




with constant term 


3.8(2.2) 


1.97(7) 


7.5(4.5) 


4.32(95) 


2.58(25) 


-5.58(42) 


-0.51(7) 



To compare chiral extrapolations with experiment, we have also plotted selected experimental data in Fig. [71 As 
noted in the introduction, there are still unresolved experimental questions, and we have indicated the range of possible 
values of {r\) 2 that can be extracted from present experiments by showing two extreme results from the literature. 



We note however, that in Ref. fUjl it was already observed that the leading one-loop SSE formula for (rj) 2 (Eq. is dominated by 

the leading chiral logarithm and dropped below the level of the lattice data available at that time for values of the pion mass as low 
as < 200 MeV. This prompted the authors of Ref. Il6l to exclude the isovector Dirac radius from the simultaneous fit. Likewise, 
the authors of Ref. [l7ll obtained huge, unrealistic values for the isovector Dirac radius when trying to enforce a fit of the logarithm- 
dominated behavior onto their data. Given these two negative precedents, we consider our "fit" to the isovector Dirac radius data to 
be of exploratory nature, testing the limits of applicability of the leading one-loop SSE results given in Eq. (142 1 . 
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The highest value is from PDG 2008 [53J and the lowest value is from a dispersion analysis including meson continuum 
contributions f(3(| ■ We note that none of the chiral fits simultaneously yields a good fit to the lattice data while also 
agreeing with experiment within statistical errors. 

To see how strongly the lattice results deviate from the SSE formulae, we also try to determine some of the low- 
energy constants using experimental results at the physical pion mass. We use the values in Tab. [V] as input, and also 
set ca = 1-5 and cy = —2.5 GeV^ 1 . Now for (r^) 2 , we have only the counter-term B[ to determine. Constraining the 
curve to go through the higher experimental value of {r\) 2 = 0.637 fm 2 gives B\ (\ = 600 MeV) = 1.085, resulting in 



the solid curve shown in Fig. 8(a) For comparison, we also plot the dashed curve that is fixed to go through the lower 
experimental value (r") 2 . The curve rises much more rapidly than the lattice data as the pion mass decreases. From 
the slope of the leading one-loop SSE curve near the physical point and the weak pion mass dependence displayed by 
our data we estimate that the applicability of Eq. (|42[) for (r\) 2 may be much less than 300 MeV. 

Without the constant term in Eq. (f5"5)) . k v ■ {r^) 2 does not have any free parameters, which yields the solid curve 



in Fig. 8(b) The curve undershoots the physical point by about 5%, which may be well accounted for by the 
uncertainties in the chosen values of the low-energy constants. Including the higher-order term C of Eq. (|55|) can of 
course shift the curve up to exactly reproduce the product of physical Pauli radius and anomalous magnetic moment. 
However, the departure of the quark-mass dependent curve from the lattice data displayed in Fig. |8(b)| indicates that 
the leading one-loop SSE formula for k v ■ (r^) 2 of Eqs. (|43l [55)) should only be trusted for pion masses much less than 
the currently available 300 MeV. Judging from the steep slopes displayed by both the curves for the Dirac and Pauli 
radii as opposed to the almost mass-independent nature of the lattice data, it is conceivable that the lea ding one-loop 
SSE formulae may only be applicable at pion masses well below 300 MeV, as already suggested in Ref. fig . 

The anomalous magnetic moment still has two free parameters, E\ and In addition to the physical point, 
we need another data point to determine both parameters. We choose to use our = 355 MeV result in the 
determination, since this point is the most accurately calculated and its relatively large pion mass makes it less 



susceptible to finite volume effects. The resulting curve (the solid line) is given in Fig. 8(c) For comparison, we also 
show the curve using the leading-order SSE formula in Eq. |39|) (the dashed line). In this case, only the experimental 
point is included to determine k° . We can see that the dashed line deviates greatly from the lattice data. This is not 
surprising, as the dominating contribution to n v is the term linear in m ff , the coefficient of which is determined by 

3 f p2 . This is clearly not the case in our data. Regarding the limit of applicability of Eq. (|40|) (which includes the 

dominant next-to-leading one-loop corrections to the strict C(e 3 ) SSE result of Eq. (f3"5)) ). the plot in Fig. 8(c) does 



not give us a clear indication up to which pion mass the formula can be quantitatively employed. Furthermore, we 
observe that the "normalized" anomalous magnetic moments display a flat pion-mass dependence around 2.5 nuclear 
magnetons. The new dynamical DWF data extend this "plateau" of the normalized magnetic moments — which was 
already observed at much larger pion masses in the quenched simulation of Ref. (l6j — now into the region of pion 
masses as low as 300 MeV. Surprisingly, we can find no indication of a rise in the magnetic moment at these low pion 
masses, although the onset of such a rise had been anticipated for pion masses around 300 MeV in the fit results of 
Ref. [II (see Fig. 11). 

Overall, these curves show much stronger curvatures than our lattice results. Even with pion masses as light as 
300 MeV, the C(e 3 ) SSE formulae do not seem to be consistent with our data. There are several possible explanations 
for the inconsistencies. One is that the pion masses in our simulations are still too heavy for the SSE formula at this 
order to be applicable, and the higher-order contributions may not be negligible in this range. The other possibility 
is that our results still suffer from uncontrolled systematic errors, such as finite volume effects, especially at the light 
pion masses. This will be discussed later in Sect. [V] We want to point out that our limited number of data points is 
not sufficient to constrain the chiral fits, which clearly demonstrates the need for calculations at lighter pion masses. 
Thus we do not regard our results in Tab. I VII as conclusive. Rather, we take it as an indication of the difficulty of 
chirally extrapolating currently available lattice data. 

Also plotted in Fig. [5] arc our domain wall results at = 330 MeV at a coarser lattice spacing [37} (a ~ 0.114 
fm), as well as our updated mixed-action calculations [27} at a lattice spacing of about 0.124 fm. These results are 
roughly consistent with the fine domain wall results, indicating that the discretization errors may be small. 



2. Chiral extrapolations using covariant baryon chiral perturbation theory 

In this section we apply a different formulation of SU (2) chiral effective field theory in the baryon sector, without 
explicit A (1232) degrees of freedom: covariant Baryon ChPT as introduced in Ref. [6l| with a modified version 
of infrared regularization (Ji?-scheme). For details about the formalism and differences from the standard infrared 
regularization introduced by Becher and Lcutwyler [34| . we refer the reader to Refs. [HI, [H, [62j. The expressions 
for the ra^-dependence of the mean squared isovector Dirac and Pauli radii and the isovector anomalous magnetic 
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Figure 7: Chiral extrapolations for the isovector Dirac radius, Pauli radius and the anomalous magnetic moment using the 
C(e 3 ) SSE formula with (solid curves) or without (dashed curves) the constant term in Eq. (|55|l . In both cases, (r\) 2 and 
K v ■ ( r 2) 2 are fit simultaneously, while k„ is fit separately with ca determined from the simultaneous fit. 



moment have been derived in (33j up to order p 4 , i.e. at the next-to-leading one- loop accuracy and are collected 
below 9 . 

For the isovector mean squared Dirac radius, the expression is given as 



K) 2 = B d + [W) 2 ] (3) + [K) 2 ] m + 0K) 



I (4) 



(56) 



In Ref. [33l | the form factor slopes an< I P2 are use( I, which are related to our notation for r" and rJJ by p\ = ^(^i) 2 and p\ = ^/^-(r^) 2 
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Figure 8: SSE chiral fits constrained to go through the physical points using the input in Tab. [V] as well as ca = 1.5 and 
cv = —2.5 GeV - . The mixed-action results at m n = 355 MeV are shifted slightly to the right for clarity. In (a) the solid 
curve is constrained to go t hrough the physical result given in PDG 2008, and the dashed curve is constrained to go through 
the result given in Ref. [6C|. In (b) the curve is drawn using the input low-energy constants according to Eq. (|43p . In (c) the 
solid curve is constrained to go through the physical point as well as our DWF result at = 355 MeV using Eq. (|40[) . while 
the dashed curve is constrained to go through the physical point using Eq. (|39|l . 
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where 



B cl = -124(A), 



I (3) 



1 



16tt 2 F 2 M 4 



7g\M 4 + 2(5g 2 A + 1)M 4 log — - + M 4 - \hg\m\M 2 

A 



g 2 A ml{lbml - 44M 2 ) log 



M 



g A m * 



1%<k 2 F 2 M 4 ^JAM 2 - m 2 

ic 6 g A ml 

l6ir 2 F 2 M$ y/AM 2 - m 2 



[15m 4 - 74m 2 M 2 + 70M 4 ] arccos [r^j) , 



m 7r (m 2 r — 3M )arccos 



2Mr 



4M 2 - m 2 n 



M 2 + (M 2 -m^log^ 



(57) 



(58) 



(59) 



The terms contributing up to and including 0(p l ) are denoted by the superscript (i). Without any loss of generality, 
the regularization scale A is set equal to Ma, the nucleon mass in the chiral limit. The low-energy constants de and Cq 
appear, respectively, in the third- and second-order ttN effective Lagrangian. The mass function M must be identified 
with M if one truncates the previous expression at 0(p 3 ), whereas at order p 4 , according to Ref. [11], M should be 
replaced by (32^ 



M A r(m 7r ) 



M - 4 Cl m; 
3mi 



ig 2 A m\ 



327r 2 F2 A /4 



Ml: 



Ml 



4ci 



M 3 



2M n 



128tt 2 F2 
-4ei(A)m 4 



65i 
M 



C2 



3gig|rn| 
8ir 2 F 2 M 2 8 V 







8ci + c 2 + 4c 3 log 



A 



(60) 



where c\, c 2 and C3 are second-order low-energy constants and e£(A) denotes an effective coupling consisting of a 
combination of fourth order low energy constants. In our current analysis, we always include terms up to 0(p 4 ), 
hence M in all the CBChPT expressions presented here should be identified with M^(m^). 
The pion mass dependence of the isovector Pauli radius is given by 



K v (m^) ■ {rlf = ^ (b c2 + [k v ■ (r v 2 ) 2 }^ + [k v ■ (r«) 2 ]( 4 >) + 0(m,), 



(61) 
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where 10 



B c2 = 24M e£ 4 (A), 
[*v WW® - 



9 2 A M 



\%v 2 F 2 M^{m 2 - U-I 2 ) 



124M 6 + 105mlM i - 18m 4 M 2 



6(3m^ - 22Af 2 m 4 + 44M 4 m 2 - 16M 6 ) log — 



9a M o 



8Tr 2 F 2 M 5 m^(4M 2 - ml) 3 / 2 



M 



9mi - MM 2 ml + 246Af 4 m 



f4 4 



[k v ■ (rj) 



^21 (4) 



216A/ 6 m 2 + 16A/ 6 



9 2 A c & m l 



arccos 



16tt 2 F 2 M^(4M 2 - mlfl 2 
1 



\2MJ 1 

[4m 4 ~ 27m l M o + 42M o] arccos 



( U| 



16n 2 F 2 M$(ml - 4M 2 ) 



16c 4 M ' + 52<^M e - 4c 4 m^Af b - 14c 65 im 2 Af 4 



13 5 2 m 2 M 4 + 8(3 5 2 - c 4 M )K ~ 4M 2 )M 4 log + 4c 6 < ? > 4 r M 2 



5 i(m 2 - 4M 2 )(4c 6 m 4 r - 3c 6 m 2 M 2 + 24M 4 ) log -f 

Mo 



For the isovector anomalous mag netic moment, the C(p 4 ) CBChPT expression is 

Mn 



M 



c 6 - 16M mle r 106 (X) + Sk^ + 6k<£> + C(mJ), 



where 



(3) _ glm^Mp 
8ir 2 F2M 3 



(3m 2 -7M 2 )log^-3M 2 



g 2 A m^Mo 



8tt 2 F 2 M 3 ^/AM 2 



[3m 4 - 13M 2 m 2 + 8M 4 ] arccos 
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32n 2 F 2 M 2 
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4M 2 (2c 6 <?i + 7fi + c 6 - 4c 4 M ) log ■ 
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2Mr 



(62) 



(63) 



(64) 



(65) 



(66) 



(67) 



Note that -^cg is equivalent to k° in Ea.(|4"D]l. 

In our chiral extrapolations, we treat gA, F w , C2, C3 and c 4 as input parameters. The available information about 
the chiral limit values of q a and F^ have been discussed in the previous section. We set the second-order couplings 
consistently with Refs. (r33l. Im. I65T] 1 



We summarize these values in Tab. I VIII 
Table VII: Input values for the covariant baryon chiral fits. 



g A F-jy [GeV] 


c 2 [GeV- 1 ] 


c 3 [GeV" 1 ] 


c 4 [GeV" 1 ] 


1.2 0.0862 


3.2 


-3.4 


3.5 



10 We note that C in Eq, l|55j l is equivalent to eL (A) 

11 For a discussion about the value of C3 see [66, 67] 
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We determine Mq, c\ and eJ(A) appearing in Mjv(m^) by fitting the nucleon masses from the three fine DWF 
ensembles to Eq. (|60p . The fit values are tabulated in Tab. IVIIII and the resulting fit curve is shown in Fig. [51 The 
fit (denoted as "Lattice only" in the table) is in excellent agreement with the physical nucleon mass, but the small 
number of data points included in the fit gives substantial statistical errors. To better constrain the value of Mq, 
which is needed in the subsequent fits, we also fit the data with the experimental point as a constraint (denoted as 
"Lattice+Exp."). The results are again shown in Tab. IVIHl The two fits give consistent results, and we will use 
central values of Mq, c\ and e\(X) determined from the "Lattice+Exp." fit subsequently 

For comparison, we also plot the coarse (a = 0.114 fm) domain wall result at w 330 MeV, as well as the mixed- 
action results [5l| at a = 0.124 fm in Fig. [9j We see that these results are qualitatively very consistent, indicating the 
discretization errors are small. 



Table VIII: Low-energy constants determined from the fit to the pion mass dependence of the nucleon mass using the C(p 4 ) 
CBChPT expression. Only the domain wall results on the fine lattices are included in the "Lattice only" fit, and in the 
"Lattice+Exp" fit we impose that the curve goes through the physical point. 



Fit 


M [GcV] ci [GeV- 1 ] el (A = 1 GeV)[GcV" d ] 


Lattice only 
Lattice + Exp. 


0.883(79) -1.01(26) 1.1(1.3) 
0.8726(29) -1.049(40) 0.90(32) 



1.6 



DWF Results, a = 0.084 fm 
DWF Results, a = 0.1 14 fm 
Mixed-Action Results, a = 0.124 fm 
Experiment 




Figure 9: Chiral extrapolation for the nucleon mass using the 0(p 4 ) CBChPT formula in Eq. ([60]). The solid line is the 
fit to only the fine domain wall data (solid circles). The square is the coarse domain wall result, and the diamonds are the 
mixed-action results from Ref. [5lj . 



We determine the remaining four low-energy constants, cq, dg(A), ey 4 (A) and e^ 06 (A), from a simultaneous fit to 



(r%) 2 and k v using 0(p 4 ) CBChPT expressions presented previously, with the results shown in Tab. IIXI 



The large x 2 /dof value indicates that the C(p 4 ) CBChPT does not describe our data either. We compare the chiral 
extrapolations using both the CBChPT formula and the C(e 3 ) SSE formula in Fig. [TUJ The solid curves with error 
bands are the results of the CBChPT simultaneous fit, and the dashed curves are the SSE fits using Eqs. (|4"2"|) , (j4"3"|) 
and (H|) as described in Sect. lIIICll It appears that both the SSE and CBChPT expressions are not compatible with 
our data, but since many of the low-energy constants in CBChPT are fixed from phenomenology or the nucleon mass, 
the fit is better constrained than that using the 0(e 3 ) SSE expressions. This is especially important for k v , for which 
the SSE expression involves more parameters than currently available lattice data. Nevertheless, both formulations 
fail to describe our data at this mass range. 



Table IX: Fit parameters for the simultaneous fit to (rl)' 1 
have set A = Mq. 



(r^) 2 and k v using the 0(p 4 ) covariant baryon formula. We 



X7dcf c 6 4(A) [GeV- a ] 


e; 4 (A) [GeV^] 


el 06 (A) [GeV" J ] 


7.3(2.4) 4.290(46) 0.839(7) 


1.350(45) 


-0.132(37) 
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Figure 10: Simultaneous fit to (t"i) 2 , k v ■ (r^) 2 and k v using the covariant baryon formula (solid lines). The dashed lines show 
the SSE fits without the constant term for k v ■ (r^) 2 - 



IV. ISOSCALAR FORM FACTORS 



Since we have not calculated the disconnected contributions to the three-point functions for the form factors, in 
this section we give results for the isoscalar form factors as defined in Eq. ([5]) from the connected diagrams only. The 
renormalized results (using the renormalization factors discussed in Sect. IIII A[) in terms of the quark flavor content 
F i,t d (Q 2 ) = F i,2(Q 2 ) + F i,2(Q 2 ) = 3 F i,2(Q 2 ) are presented in Tab. IXlTTl IXVl First, we study the Q 2 dependence of 
both the isoscalar Dirac and Pauli form factors using phenomenological models, and then discuss briefly the chiral 
extrapolations of the results. 



A. Q dependence 



Unlike the isovector Dirac form factor. 



F i +d (0) is not set to the known value of 3. Thus we perform dipole fits to 
F^ +d (Q 2 ) separately to each ensemble using the formula in Eq. (|30|) . Similar to the isovector case (see Sect. IIII Bl) . 
the dipole Ansatz describes the data reasonably well at small Q 2 values, typically below 0.6 GeV 2 . As large Q 2 values 
are included in the fit, the fit quality becomes worse, but the fit parameters do not change significantly. Furthermore, 
the fitted values of F 1 tl+d (0) are very consistent with the expected value of 3. 

To demonstrate the quality of the fits, in Fig. [11] we show the dipole fits to all the Q 2 values. One can see that 
the data are reasonably well described by the fit curves. Also plotted is the phenomenological fit to experimental 
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data using the parameterization in Ref. (48j . although we note that no error estimate is provided and the empirical 
analysis involves many potential systematic errors discussed in the introduction. To determine the isoscalar mean 
squared Dirac radii, we follow the same reasoning as in Sect. IIII Bl and obtain them from the dipole fits with a cut at 
Q 2 < 0.5 GeV 2 . The results are shown in Tab.0 

In experiments, the isoscalar Pauli form factor shows a notable bump at Q 2 ss 0.4 GeV 2 (solid curve in Fig. [P2]) . 
although again there are no error estimates. Our data are too noisy to distinguish this feature at this moment. In 
fact, the results, shown in Fig.[T2J are rather flat. We show the constant fits to each ensemble separately, and find that 
the constants are consistent with zero within two standard deviations. The error band corresponds to the constant 
fit to the m v = 297 McV data. 

If we restrict the fits to only the small Q 2 region (< 0.5 GeV 2 ), we are able to perform linear fits to the data and 
obtain both k s ■ (r|) 2 (from the slope) and k s (from the intercept), the results of which are also shown in Tab. IXr 2 . 

Table X: Results for the isoscalar Dirac and Pauli mean squared radii and the anomalous magnetic moment. A dipole fit with 
a Q 2 cutoff at 0.5 GeV 2 is used to determine (rf) 2 . Linear fits to F£ with Q 2 < 0.5 GeV 2 are used to determine k s ■ (r|) 2 and 
« s . The results shown below have been normalized to the physical nuclear magneton. 



m n [MeV] 


xVdof (rlY [GeV" 3 *] 


XVdof 


norm / s\'Z T/~i -v r — 'Zi 

k s ■ (r 2 ) [GeV ] 


norm 


297 


0.12(35) 11.00(13) 


3.3(2.1) 
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Figure 11: The isoscalar Dirac form factor, F^ +d (Q 2 ), with dipole fits. The thick solid (red) curve is a phenomenological fit to 
experimental data [4Sl |. 



B. Chiral extrapolations 

1. Chiral extrapolations using 0(e 3 ) small scale expansion 

As is well known in ChPT (e.g. see the discussion in chiral dynamics in the isoscalar form factors of the 

nucleon starts at the 3-pion cut, i.e. at 2-loop level, corresponding to C(e 5 ) in the power-counting of SSE. Hence, 
although the C(e 3 ) SSE expressions for the pion mass and momentum transfer dependence of the isoscalar Dirac and 



Like in the isovector case, the anomalous magnetic moment quoted here is normalized to the physical nuclear magneton according to 
Eq.ltill. 
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m = 297 MeV 
= 355 MeV 
m i = 403 MeV 
Phenomenology 




Q 2 [GeV 2 ] 

Figure 12: The isoscalar Pauli form factor, _F^ +ti (Q 2 ), with constant fits. Only the error band for the fit to the = 297 MeV 
ensemble is shown. The thick solid (red) curve is a phenomenological fit to experimental data [48| . 



Pauli form factors have also been derived in [4l| and given as 



1 + B 1 



Q 2 



(47TF,) 2 ' 



F^Q 2 ) = k s , 



they cannot be utilized for chiral extrapolations. Therefore, in this section, we simply extrapolate linearly in m 2 , 
the mean squared Dirac radius to the physical point. This is shown in Fig. 1131 where we can see that the linear 
extrapolation gives a result at the physical pion mass which is much lower than the empirical value. Similarly, we 
perform a linear extrapolation for k s ■ (r^) 2 , which is shown in Fig. 1141 
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Figure 13: The isoscalar Dirac radius and linear extrapolation. The star indicates the phenomenological value obtained in 
Ref. H. 

For k s beyond order e 3 , additional terms arise including a term linear in the quark mass. Following Ref. [49| . we 
write 

k s = k° s - 8E 2 M N ml, (68) 
where k° and E2 are two unknown LECs. This linear dependence describes our data well, as is shown in Fig. 1151 
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Figure 14: The isoscalar Pauli radius and linear extrapolation. The star indicates the phenomenological value 18.4 GeV 
obtained in Ref. [f|. 
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Figure 15: The isoscalar anomalous magnetic moment, k s and linear extrapolation. The star indicates the experimental 
value [53l ). 

2. Chiral extrapolations in covariant baryon chiral perturbation theory 

The CBChPT formulae up to C(p 4 ) for (rf) 2 , (r s 2 ) 2 and k s have also been derived in [33[. We collect them here 
for completeness. We note, however, that the next-to-leading one-loop CBChPT results for the isoscalar form factors 
of the nucleon as presented in this section — just as in the case of the leading one-loop SSE-analysis discussed in the 
previous section do not contain their dominant chiral dynamics arising from the 3-pion cut. Such effects would only 
become visible at the two-loop level, i.e. starting at 0(p 5 ) in CBChPT. The results presented here are therefore to 
be interpreted with care, as several important contributions with potentially large impact on the chiral extrapolation 
functions are not included at this order. For the isoscalar mean squared Dirac radius, the CBChPT expression is 
given by 



(^) 2 = S c s 1 +[(r i ) 2 ] (3) + [(^) 



I (4) 



(69) 
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where 
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Here, again, when the expression is truncated at 0(p 3 ), M should be identified with Mo, while at 0(p 4 ), it should be 
replaced by M^{m v ) in Eq. (|6"0)) . Similarly, for n s ■ (r|) 2 , we have 

Mn 
M 

with 
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The CBChPT expression for the isoscalar anomalous magnetic moment is written as 
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As in the isovector case, we use the values in Tab. I Vill as input in the extrapolations, leaving d?, and e^ 05 (A) 
as free parameters. Since (rf) 2 , k s ■ (r^) 2 and k s all contain the low-energy constant naively we should perform 
a simultaneous fit to all three quantities, as we have done for the isovector case. However, as stated earlier, the 
dominant chiral dynamics for the isoscalar quantities appears at 0(p 5 ). We do not expect these 0{p 3 ) expressions 
to describe our data. In fact, the simultaneous fit to these three quantities gives a % 2 /dof of about 9 (see Tab. IXI|) . 
showing the difficulty in fitting these quantities consistently. Looking closely at each quantity separately, we find that 



independent fits to (rf 



,s\2 



„s\2 



and k s lead to an inconsistency in the estimation of the common parameter k", 



as shown in Tab. IXIII For demonstrative purposes, we compare the resulting fit curves from the simultaneous fit and 
the independent fits in Fig. 1161 from which we see that the independent fits provide reasonable extrapolations for the 
data, while the simultaneous fit misses the data points badly, indicating inconsistencies of the CBChPT expressions 
at this order. We also note that the extrapolated value for (rf ) 2 at the physical pion mass is about 20% lower than 
the phcnomcnological value. These observations lead us to conclude that the CBChPT expressions at 0(p 3 ) are not 
applicable in the pion mass range of our calculation. Of course, since we have not included the disconnected diagrams 
in our calculations, there are uncontrolled systematic errors which may affect the pion mass dependence. Further 
investigations are required to draw definitive conclusions for these isoscalar quantities. 
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Table XI: Fit parameters from the simultaneous fit to (V^) 2 , k s ■ (r|) 2 and k s using Eqs. (|69|) . (|720 and (|76[) . 

xV dof d 7 e 54 ei 05 (A = M ) 

8.5(2.6) -0.172(23) -0.458(24) -0.0159(41) 0.598(26) 



Table XII: Fit parameters from independent fits to (rl) 2 , k s ■ (r|) 2 and k s using Eqs. (|69|l . ([720 and (|76[) . 
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V. SYSTEMATIC ERRORS 

A. Effect of the excited states 

The correlation functions may have systematic bias due to the excited and/or unphysical oscillating states [2~il . 
|28[ [68[. To control it, we solve the overdetermined system separately for each location of the operator and examine 
the plateau for the form factors. Examples are shown in Fig. 1171 Due to the tuning of the quark sources, the 
contaminations from states other than ground are suppressed and do not contribute to the matrix element plateaus 
close to their centers. 

To put quantitative bounds on possible bias, we study the excited states in the nucleon correlators. The nucleon 
two-point correlation functions have very precise information on the presence of the non-ground state contamination. 
For example, with our current statistics the parameters of a fit with three states are well constrained: 

C 2pt (t;P) = Z Q (P)e- Eot + Z x {P)e~ E ^ + (-l)% sc (P) e - £ °-<, Z 0A > 0. (79) 

Having estimated the energy gap AEiq(P) = E\(P) — Eq(P) and the magnitude of the contamination Z\(P)/Zo(P), 
one can put bounds on the excited state contribution to the matrix elements computed from the two- and three-point 
lattice nucleon correlators. 

The ratio formula (|23| for physical matrix elements has two factors: R v = RnRa- Excited states can potentially 
contribute to cither one. First, we study the asymmetry ratio, Ra, defined in Eq. (|22p . As was pointed out above, 
this factor compensates the asymmetric r-dependence in i?jv, and in the absence of excited states it would be equal to 
exp \—{E' — E)(t — T/2)]. Although this factor involves different two-point functions, their excited state contributions 
appear to cancel each other to a large extent, as shown in Fig. [18] The left panel of Fig. [18] shows the ratio of Ra to 
the exponential result in the absence of excited states 



„ , . C 2p t(T-T,P)C 2pt (T,P>) 

Ra(t) _ V C 2p t(T-r,P')C2 P t(r-P) / go s 



e -(E'-E)(r-T/2) e -(E'-E)(r-T/2) 

where (E' — E) in the denominator is determined by the best fit to Ra in the range 3 < r < 6 . The fact that this 
ratio is unity to within 1% over a plateau from 3 < r < 9 indicates that excited state contributions are negligible. 
Furthermore, the right panel of Fig. 1181 shows the effective ground state energy difference 



5E cit {t) = log 



C 2p t(i,P') , C 2pt (t,P) 



C 2pt (i + l,P') C 2pt {t + l,P) 



(81) 



which in the absence of any excited state contaminants, would simply be SE cS (t) = (£" — E). For comparison, the 
values of E 1 — E determined above are plotted on the same graph, and agree nicely in the fiducial range 2 < r < 10. 
Thus, we neglect small contaminations from this factor. 

Second, we estimate the contribution to Rn defined in Eq. (|2ip assuming only one excited state and no oscillating 
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Figure 16: 0(p 4 ) CBChPT fits to (rf) 2 , n a 
resulting fit parameters summarized in Tab. 
Tab. EH 



IxTT 



2 and k s . The solid lines are fits to the three quantities separately with the 
The dashed lines are simultaneous fits with the parameters summarized in 



term 13 : 

C 3p t(r,T)w C 3pt (r,T)\ 

C 3p t(r,T) 
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l + ^i2io(r) + ^i2i (T-r) 
+ ^^loCrJKoCr - r) - \ (SR n + 6R' n ) 



(82) 



7 {>) 

^1 -A£ (/) T 
Z 



8R§(T/2) 



(83) 



We neglect the contribution of oscillating states because they decay even faster than excited states. 
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Figure 17: Nucleon form factor plateaus for the lightest = 297 MeV ensemble. 
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Figure 18: The left panel shows the ratio A(t) in Eq. (|80[1 . The right panel shows the effective energy difference (|81l) in lattice 
units and the fit values of E' — E used in the left figure. The degree to which the contaminations to all P' ^ two-point 
correlators are canceled by the contamination to the P = correlator is remarkable 



and we have expanded to leading order assuming that 8R{{ < 1. The value of the suppression factor 6R§(t) IS 
shown in Fig. 1191 Its values arc estimated from the fit parameters Zq.i, -Ea.o in Eq. (|79[) . and the errors are computed 

using the Jackknife procedure. Note that 5i?^g(r) falls off steeply with r. As a result, its contribution can be easily 
detected and removed by fitting the plateau with 

R°(r) ^C + C ie - AE + C[e~ AE '( T -^ (84) 

From Fig. [15] one may estimate the last two terms in the contamination formula (fB"2)) . suppressed by 6R±l and 
SRxo(t)SR' w (T — t). If one further assumes that the excited state matrix elements are at most of the same order as 
the ground state elements, < 1, the effect of the last two terms in Eq. ([82]) is well below 1%. It is also worth 
noting that higher momentum matrix elements with p = (0, 0, 2) would contain substantially larger contamination, 
as compared to lower momenta. Such matrix elements are excluded from our analysis. 

Finally, we compare the form factors extracted using the plateau average and fitting the r-dcpendence to Eq. (|84[) . 
Due to the uncertainty in the two-point correlator fitting parameters, we perform fits for a range of mass gaps 
AMjy = 0.4, 0.6 and 0.8, which bracket the fitted values from different fitting ranges and fitting with or without 
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Figure 19: Suppression factor for the excited state contributions SR U )(t) (|83|) . as estimated from fittingthe two-point function. 
Note that the actual matrix element do is not included in the plotted value, which therefore shows only the relative fall-off 
of the exponential tail contamination. Note also that the factor for the p = (0, 0, 2) state is substantially larger than for the 
other states shown. 

the oscillating term in Eq. ([75)) . The energy gaps AE for the P / states are computed using the continuum 
dispersion formula. The result is statistically independent of the mass gap value used (see Fig. [20]) and is stable when 
fitting inside the region 2 < r < 10. The complete consistency between conventional plateau averages and results for 
which excited state contaminants are explicitly included in the analysis and separated from the physical ground state 
contribution clearly indicates the absence systematic errors from excited state contaminants in our present results. 

In addition, we have also compared results with two different source-sink separations, T = 12 and T = 14. If the 
coherent sink technique were ever to introduce additional noise into the calculation, one would expect it to be worst 
for the larger T, for which the first adjacent unwanted sink is closer. Hence, in the case of T = 14, we have used 
independent sinks to check that this is not a problem. One of the typical plateaus comparing T = 12 and T = 14 
separations is plotted in Fig. [21] and shows agreement within statistics. Separations 12 and 14 are also compared in 
Fig. [20] where we show each of the form factors computed on a subset of the am q = 0.004 ensemble using independent 
backward propagators and the larger source-sink separation T = 14. The agreement of results that use two different 
separations and techniques directly indicates that our method does not suffer from the systematic effects due to 
excited states or the coherent propagator technique. 

B. Finite volume dependence 

Ideally, we would like to control systematic errors arising from volume dependence using an effective field theory 
that describes the dependence of the observables of interest as a function of spatial volume and pion mass, and a 
set of calculations of the lattice observables with a specified action for a range of volumes at pion masses for which 
the effective field theory is applicable. Verifying that the effective theory fits the measured volume dependence 
with low energy constants that are consistent with other lattice and phcnomenological constraints would then assure 
solid theoretical and computational control of finite volume effects. To date, this program has not been carried out 
completely for form factors 14 with any lattice action, and finite volume effects have been a convenient excuse for any 
disagreements with experiment. Hence, it is useful to examine the available data for domain wall fcrmions and to 
assess what quantitative evidence there is for or against significant finite volume corrections. 

To examine volume dependence carefully, it is important to only compare lattice calculations at different volumes 
that use precisely the same action and computational methodology. In this context, we believe it can be seriously 
misleading to argue on the basis of plots containing a variety of calculations with different actions, analysis techniques, 
renormalization schemes, etc. In addition, we will find it useful to distinguish between forward matrix elements, for 



We note that the Regensburg group has recently started the extension of the C(p 4 ) CBChPT calculation for the isovector form factors 
of the nucleon to a finite volume in the p-regime l69t . 
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Figure 20: Comparison of form factors extracted from plateau averages and from fitting plateaus with formula (|84|l for the 
ensemble with the lightest pion mass = 297 MeV. The result is stable with variation of the mass gap AA/jv, which means 
that the contamination is small. All but the last group of points use source-sink separation T = 12. The last group, calculated 
for 330 gauge configurations, uses the larger separation T = 14 and independent backward propagators. Each form factor value 
is divided by the central value of the dipole fit. Tab. [II] lists the momentum combinations corresponding to each index on the 
horizontal axis. 
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Figure 21: Comparison of F™ d plateau using coherent backward propagators with T = 12 and independent backward propa- 
gators with T = 14. The momentum transfer Q 2 corresponds to (000|011). 



which we have applied detailed effective field theory formulae to finite volume corrections [55j, and off-diagonal matrix 
elements, for which we have not yet done so. Due to the high computational cost of domain wall fcrmions, the only 
high precision calculations of nucleon observables we are aware of with two large volumes for light quarks are the 
mixed action calculations for 355 MeV pion mass in volumes of spatial extent 2.5 fm and 3.5 fm, corresponding to 
ra^L = 4.4 and 6.2 respectively (2(|[27]]- Thus, we will base our computational arguments on these results. We begin 
with the forward matrix element corresponding to the axial charge, which was first calculated for these two lattice 
volumes in Ref. (55[. In that work, the values of the axial charge calculated at the two volumes agreed within their 
statistical errors of 5%, which when combined quadratically indicated that the volume dependence was less than 7%. 
A new high statistics calculation (2?| yields the calculated fractional difference between gA at 3.5 fm and 2.5 fm of 
(—0.4 ± 1.9)% at ni-x = 355 MeV. In effective field theory, the loop integrals over momenta arising in an infinite 
volume are replaced by sums over discrete momenta in a periodic finite volume, and explicit expressions are available 
for the axial charge in a finite volume (57l l70l . l7ll ]. Using the low energy parameters determined from fitting the 
data in Ref. (55[, chiral perturbation theory specifies that the fractional difference between gA at 3.5 fm and 2.5 fm is 
+0.4%. This is consistent with the observation in Ref. [24[ that for a reasonable value of ca ~ 1-5, the finite volume 
corrections are quite small. Hence, for the forward matrix element gA, all the evidence is consistent in indicating that 
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the systematic error arising from a spatial box of extent 2.5 fm is of the order of one percent. 

Since we have not performed a similar analysis of chiral perturbation in a finite volume for form factors at finite 
momentum transfer, our only recourse at present is comparison of the mixed action numerical results. Because we 
are most interested in chiral extrapolation of rms radii, we have focussed on the most accurate means of calculating 
the slope at zero momentum transfer. Thus, we take the momentum combination (1, 0, 0)2n/L, (0,0,0), on each 
lattice and fit the form factor at the resulting momentum transfer with a one-parameter dipole formula, from which 
we determine the slope at the origin. Comparing the results for r\ and quadratically combining the errors for the 
two independent calculations, we find that the fractional difference between r\ at 3.5 fm and 2.5 fm is (3.7 ± 2.6)% at 
TOjt = 355 MeV 15 . Further evidence suggesting finite volume corrections to radii are small is the fact that in Ref. [29| . 
even decreasing the lattice size to 1.8 fm yields small changes in r\ and r%- Both by virtue of the fact that Monte 
Carlo calculations of the slope are intrinsically noisier than for forward matrix elements and fact that we have not 
performed an effective field theory analysis of the finite volume corrections, our control of the volume dependence of 
form factor radii is worse than for gA- Whereas the error may also be on the order of one percent as in gA, we cannot 
completely exclude a result at the upper limit of the error bars of the order of 6%. Thus, presently, it cannot be 
excluded that a rapidly growing finite-size effect below = 355 MeV is affecting the pion mass dependence of r\ in 
our data. This will have to be resolved in the future. 



VI. COMPARISON WITH PREVIOUS CALCULATIONS 



We briefly compare our results and conclusions with those of previous calculations. We start with the isovector 
Dirac and Pauli radii, r\ and r". Previous calculations using Wilson fermions had reached pion masses of about 
400 MeV. Both the quenched [l|| E3 an d Nf = 2 unquenched results showed a mild pion mass dependence 
for r\ and r\. A recent calculation on Nf = 2 + 1 domain- wall fermion configurations at a = 0.114fm extended 
the range of pion masses down to 330 MeV [2j|. These results show that the very mild upward trend of r\ and r% 
extends down to that pion mass. Summary plots comparing results for Wilson fermions with zero and two flavors 
and domain- wall fermions for zero, two, and 2+1 flavors are given in Figs. 14 and 19 of Ref. (2a |. For the case of r\, 
which has smaller statistical errors, for each action, the data tend to lie on straight lines with comparable small slope 
and some scatter in normalization, with perhaps a hint that the Nf = 2 calculations, performed on box sizes 1.9 fm, 
lie somewhat low. The r% data also appear to lie on straight lines with similar small slope, albeit with larger scatter. 
Our results for r\ and r^, which extend down to = 300 MeV, also show a small pion mass dependence, and are 
consistent within statistical errors with the 2+1 flavor domain wall results on a = 0.114 fm lattices. We conclude 
that this flat behavior, surprising as it is from the chiral effective theory point of view, is genuine. The one-loop SSE 
formulae of Sect. IIII C II cannot accommodate this "flat" pion mass dependence in the radii down to such low values 
of the pion mass ~ 300 MeV, with or without the inclusion of a higher order "core" term. Indeed, the curves shown 
in Figs. 8(a)[ |8(b)| indicate that the SSE calculation would have favored an upward trend in the extracted isovector 



radii which should have become visible in the pion-mass range studied in this work, consistent with the expectations 
drawn in Ref. [l|. The only explanation for this behavior available at the moment is that the leading one-loop SSE 
calculation is only valid for pion masses < 300 MeV. 

As for the anomalous magnetic moment k v , our results are in very good agreement with those obtained with Nf = 2 
dynamical Wilson fermions in [2l|, with recent Nf = 2 twisted- mass results [23| and with the recent Nf = 2 + 1 
domain- wall calculation [29[ . We remark that our Fig. [5] displays the "normalized" anomalous magnetic moment 
K norm , while Fig. 17 of Ref. [29[ shows the magnetic moment normalized by the quark-mass dependent nucleoli mass, 
K iat rpj ic c [jfj' ercnce between the two figures 16 reflects the dependence of the nucleon mass, which is quite strong 
(see Fig. All in all, for k v too, the calculated pion mass dependence is rather mild, and results at lower pion 
masses will have to bend upwards rather sharply if they are to agree with the experimental value. 



15 We note that as discussed in Ref. [3tJ , dipole fits of all the form factor data out to some fixed cutoff yield discrepancies in the dipole fits 
between the two volumes that increase as the cutoff increases, but this comparison focusses on other features of the form factor besides 
the radius that we seek to chirally extrapolate. 

16 The numerical difference between the lattice data and the experimental value is smaller in the case of fq at . 
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VII. SUMMARY AND CONCLUSIONS 

We have presented lattice calculations of nucleon form factors with Nf = 2+1 flavors of dynamical domain wall 
fermions on fine 32 3 x 64 lattices with a = 0.084 fm at pion masses of 297, 355, and 403 MeV that achieve a new level 
of precision in both statistical and systematic errors. Statistical errors have been reduced by using from 3600 to 7064 
measurements of operators at a given mass by performing 8 measurements per lattice and verifying their statistical 
independence. Statistical errors and error correlations have been carefully analyzed in our overdetermined analysis, 
which combines as many stochastically distinct measurements of the same physical form factors as practical. 

Because of the high level of statistical precision, we have carefully investigated and controlled potential sources of 
systematic error. We have ruled out systematic errors arising from the source-sink separation in two different ways. 
First, we have derived analytic expressions for the contamination by excited states, and, using lattice data from 
two-point correlation functions, have shown quantitatively that the coefficients of excited state admixtures in these 
expressions yield negligible contributions to the observables of interest. In addition, we compared explicit calculations 
with source-sink separations T — 12 and T = 14 and have shown that results from the two source-sink separations 
are indeed statistically consistent, as expected from the excited state analysis. We have verified that even in the 
worst case — the lightest pion mass and maximum source-sink separation — results calculated using the time-saving 
coherent sink technique are consistent with results calculated with conventional independent sinks. By comparison 
with companion calculations on a coarse lattice with a = 0.114 fm, we have verified that lattice spacing errors are 
small. Finally, based on the overall consistency between the recent high-statistics mixed action results and the current 
work, we have presented two arguments that finite volume corrections to the present calculations in a volume of spatial 
extent 2.5 fm are small. First, the forward matrix element gA changes by a very small amount when the spatial extent 
is changed from 2.5 to 3.5 fm: +1.0% in chiral perturbation theory and (—0.4+ 1.9)% in explicit lattice calculations. 
For the more complicated case of off-diagonal matrix elements, the measured fractional change in r\ is (3.7 ± 2.6)%. 

The high precision of the calculated form factors is shown in Figs. 14151 where, in order to see the discrepancies with 
dipole fits, we plotted the ratio of the lattice calculations to the best dipole fits on an expanded scale. This precise 
data enabled us to extract the Dirac radius, r\, Pauli radius r%, and anomalous magnetic moment k v with much 
smaller errors than in earlier calculations and to study chiral extrapolations to correspondingly higher precision. In 
contrast to earlier studies in which the lattice error bars were sufficiently large that the data appeared to be consistent 
with chiral perturbation theory, in this work we have shown that the ra, dependence of the lattice results for (r^) 2 , 
(r?;) , and k v at the three masses 297, 355, and 403 MeV cannot be simultaneously fit by either C(e 3 ) SSE or NNLO 
CBChPT. The data points for {r\) 2 rise too slowly with decreasing m n and the data for (r^) 2 are too fiat to be fit 
by either the SSE or CBChPT curves that rise smoothly with decreasing m n to approach the experimental results. 
Since there happen to be three free parameters in SSE to fit the 3 measured values of k v , the SSE can actually 
fit the anomalous magnetic moment, but CBChPT, which is physically constrained to rise with decreasing m n , is 
also seriously in conflict with the lattice measurements of n v at the accessible masses. Similarly, we were unable to 
simultaneously fit the isoscalar quantities (rf 2 ) and k s , which, to this order of ChPT, have fewer parameters. 

With the present data at these three pion masses, we see three possible explanations for the discrepancy with 
chiral perturbation theory. One possibility is an outright error somewhere in the lattice calculations. However, by 
virtue of meticulous checks, key calculations with independent codes, and the qualitative similarity of our results 
to those of other groups [ij], HH [H, HI, H9[ , we believe this is unlikely. A second possibility is that finite volume 
effects arc significantly larger than the estimates we obtained from our 355 MeV mixed action studies for spatial sizes 
2.5 fm and 3.5 fm. This possibility clearly warrants further study of chiral perturbation theory for off-forward matrix 
elements in a finite volume and careful high statistics studies in a series of volumes. The third possibility is that chiral 
perturbation theory at the present order is not applicable for this range of m^. Indeed, significant problems have 
previously been encountered in describing the dependence of baryon masses, and one observes, for example, that 
the highly linear dependence of the nucleon mass on seen in a variety of lattice calculations with different actions 
can only arise from an apparently unnatural cancellation of analytic and non-analytic terms in chiral perturbation 
theory [5l|] . This possibility clearly warrants lattice calculations at a series of lower values of m n all the way down to 
the physical pion mass. 

The last two possibilities each raise very interesting and important questions in hadron structure. Given the high 
computational cost of chiral fermions relative to improved Wilson fermions and the fact that there are no crucial 
operator mixing problems in form factors necessitating exact chiral symmetry on the lattice, it appears that the most 
expeditious means of understanding the volume dependence and behavior down to the physical pion mass will be 
with an appropriate form of an improved isotropic Wilson action. Such calculations are clearly essential for further 
progress in understanding the fundamental structure of the nucleon. 
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Appendix A: TABLES 



Table XIII: Renormalized results for the Dirac and Pauli form factors from the ami = 0.004 ensemble with m-n ~ 297 MeV. 



(aQ) 2 Q 2 [GeV 2 ] F? 



Fl 



u + d 



F u 



=7T 



Fn 



F" 



F, 



u-\-d 



F 2 " 



=7T 



0.000000 


0.000 


2.004(5) 


0.037025 


0.203 


1.457(7) 


0.037235 


0.204 


1.462(13) 


0.071421 


0.392 


1.132(9) 


0.072155 


0.396 


1.130(21) 


0.077106 


0.423 


1.089(11) 


0.103678 


0.569 


0.912(13) 


0.114341 


0.627 


0.870(13) 


0.154213 


0.846 


0.696(15) 


0.191447 


1.050 


0.591(13) 



1.004(3) 

0.683(4) 

0.681(8) 

0.497(5) 

0.491(12) 

0.470(6) 

0.376(7) 

0.350(7) 

0.256(8) 

0.204(7) 



3.008(7) 1. 
2.140(9) 0. 
2.142(17) 0. 
1.629(12) 0, 
1.621(29) 0, 
1.560(15) 0, 
1.288(18) 0. 
1.219(17) 0, 
0.951(20) 
0.795(17) 0, 



000(5) 

774(7) 1 

781(13) 

635(8) 

639(19) 

619(10) 

536(11) 

520(12) 

440(13) 

387(11) 



050(46) 

956(106) 

853(35) 

995(104) 

745(42) 

685(44) 

634(41) 

479(36) 

408(25) 



092(27) 
056(60) 
822(23) 
.693(55) 
.806(26) 
647(26) 
601(26) 
.463(23) 
388(16) 



-0.042(59) 
-0.100(139) 
0.031(50) 
0.302(127) 
-0.061(54) 
0.038(61) 
0.034(54) 
0.016(49) 
0.020(34) 



2.142(46) 

2.013(102) 

1.675(32) 

1.688(108) 

1.550(43) 

1.332(39) 

1.235(43) 

0.943(36) 

0.796(25) 



Table XIV: Renormalized results for Dirac and Pauli form factors from the ami = 0.006 ensemble with ~ 355 MeV. 



{aQ) 2 Q 2 [GeV 2 



F u 



F'i 



0.000000 


0.000 


2.000(3) 


0.037176 


0.204 


1.478(4) 


0.037348 


0.205 


1.471(8) 


0.071948 


0.395 


1.162(6) 


0.072559 


0.398 


1.155(12) 


0.077106 


0.423 


1.134(8) 


0.104730 


0.574 


0.946(9) 


0.114455 


0.628 


0.904(9) 


0.154213 


0.846 


0.736(10) 


0.191561 


1.051 


0.619(9) 



1.000(2 
0.691(2 
0.691(4 
0.509(3 
0.503(7 
0.491(4 
0.389(5 
0.365(4 
0.285(5 
0.225(4 



Fl 



u-\-d 



2.999(4) 

2.169(6) 

2.162(10) 

1.671(8) 

1.659(17) 

1.625(10) 

1.336(12) 

1.269(11) 

1.021(14) 

0.844(11) 



F u- 



F 



u-\-d 



pu-a 



1.000(3) 

0.788(4) 

0.779(8) 

0.653(5) 

0.652(11) 

0.643(7) 

0.557(7) 

0.540(8) 

0.452(9) 

0.393(7) 



1.132(31) 
1.226(71) 
0.909(25) 
0.866(62) 
0.878(28) 
0.745(28) 
0.693(29) 
0.560(25) 
0.476(17) 



-1.192(17) 
-1.130(38) 
-0.923(15) 
-0.875(33) 
-0.870(16) 
-0.727(17) 
-0.656(16) 
-0.505(15) 
-0.407(11) 



-0.061(40) 
0.096(93) 
-0.014(34) 
-0.009(79) 
0.008(36) 
0.018(37) 
0.036(37) 
0.056(33) 
0.070(23) 



2.324(30) 
2.356(66) 
1.832(24) 
1.741(59) 
1.748(28) 
1.472(27) 
1.349(28) 
1.065(25) 
0.883(17) 



Table XV: Renormalized results for Dirac and Pauli form factors from the ami = 0.008 ensemble with m^ w 403 MeV. 



(aQ) 2 Q 2 [GeV 2 ] F? 



F" 



Fi' 



Fl' 



FS 



F 



0.000000 


0.000 


2.006(3) 


0.037277 


0.204 


1.502(4) 


0.037427 


0.205 


1.499(8) 


0.072306 


0.397 


1.180(6) 


0.072839 


0.400 


1.168(12) 


0.077106 


0.423 


1.160(8) 


0.105450 


0.578 


0.965(8) 


0.114533 


0.628 


0.924(8) 


0.154213 


0.846 


0.771(11) 


0.191639 


1.051 


0.633(9) 



1.006(1 
0.706(2 
0.706(4 
0.521(4 
0.519(6 
0.508(4 
0.401(5 
0.381(5 
0.302(6 
0.234(5 



3.012(3) 

2.208(5) 

2.204(11) 

1.702(8) 

1.687(16) 

1.667(11) 

1.366(11) 

1.305(11) 

1.073(14) 

0.867(12) 



1.000(2) 

0.796(4) 

0.793(7) 

0.659(5) 

0.649(10) 

0.652(7) 

0.564(7) 

0.543(8) 

0.469(9) 

0.399(8) 



1.210(32) 
1.342(65) 
0.965(26) 
0.985(60) 
0.918(29) 
0.783(29) 
0.710(27) 
0.567(24) 
0.459(17) 



-1.193(19) 
-1.131(39) 
-0.926(17) 
-0.922(36) 
-0.861(19) 
-0.749(18) 
-0.669(18) 
-0.539(16) 
-0.442(12) 



0.016(43) 
0.211(85) 
0.038(36) 
0.063(80) 
0.058(39) 
0.035(39) 
0.041(38) 
0.028(33) 
0.017(24) 



2.403(31) 
2.473(65) 
1.891(26) 
1.908(59) 
1.779(29) 
1.532(28) 
1.379(26) 
1.106(24) 
0.901(17) 
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Table XVI: Comparison of fit Ansatze to the isovector Dirac form factors d for all three ensembles with different Q 2 cutoffs. 



ami = 0.004 





Dipole 


Tripole 


Q 2 cutoff [GeV^] 


X Vdof 


M D 2 [GeV _a ] 


x 7dof 


M~ 2 [GeV - *] 


0.3 


0.2(6) 


0.670(22) 


0.2(6) 


0.436(14) 


0.4 


0.3(6) 


0.659(19) 


0.8(9) 


0.424(12) 


0.5 


0.5(6) 


0.653(17) 


1.0(9) 


0.418(11) 


0.6 


0.4(5) 


0.652(17) 


1.0(8) 


0.417(11) 


0.7 


0.5(5) 


0.649(17) 


1.2(8) 


0.414(11) 


0.9 


0.9(7) 


0.638(16) 


1.9(1.0) 


0.404(10) 


1.1 


1.4(8) 


0.632(16) 


3.0(1.1) 


0.398(10) 



ami = 0.006 





Dipole 


Tripole 


Q 2 cutoff [GeV 2 ] 


XVdof 


ME* [GeV- a ] 


X^/dof 


M-' z [GeV~ 2 ] 


0.3 


0.5(1.0) 


0.625(13) 


0.5(1.0) 


0.407(8) 


0.4 


1.8(1.3) 


0.610(12) 


3.3(1.8) 


0.393(7) 


0.5 


2.8(1.5) 


0.602(11) 


4.8(1.9) 


0.386(7) 


0.6 


2.3(1.2) 


0.602(11) 


4.2(1.7) 


0.386(7) 


0.7 


2.1(1.1) 


0.601(11) 


3.8(1.5) 


0.385(7) 


0.9 


2.0(1.0) 


0.597(11) 


4.1(1.4) 


0.379(7) 


1.1 


2.0(9) 


0.595(11) 


4.8(1.5) 


0.375(7) 



ami = 0.008 





Dipole 


Tripole 


Q 2 cutoff [GeV 2 ] 


X 7dof 


M D 2 [GeV - *] 


XVdof 


M~ 2 [GeV~ 2 ] 


0.3 


0.09(42) 


0.592(13) 


0.09(42) 


0.386(8) 


0.4 


0.3(5) 


0.588(12) 


1.0(1.0) 


0.380(7) 


0.5 


0.9(9) 


0.582(11) 


1.9(1.2) 


0.374(7) 


0.6 


1.0(8) 


0.579(11) 


2.2(1.2) 


0.371(7) 


0.7 


0.9(7) 


0.579(11) 


2.0(1.1) 


0.370(7) 


0.9 


1.1(7) 


0.575(10) 


2.7(1.2) 


0.366(6) 


1.1 


1.0(7) 


0.575(10) 


2.7(1.1) 


0.365(6) 



Table XVII: Comparison of fit Ansatze to the isovector Pauli form factors F% d for all three ensembles with different Q 2 cutoffs. 



ami = 0.004 





Dipole 


Tripole 


Q 2 cutoff [GeV^] 


XVdof 


F 2 (0) 


[GeV- a ] 


XVdof 


F 2 (0) 


M~ 2 [GeV- 2 ] 


0.5 


1.2(1.3) 


2.89(12) 


0.820(70) 


1.2(1.3) 


2.85(11) 


0.505(40) 


0.6 


1.1(1.1) 


2.92(11) 


0.846(63) 


1.0(1.0) 


2.87(10) 


0.516(36) 


0.7 


0.9(8) 


2.93(11) 


0.847(60) 


0.8(8) 


2.87(10) 


0.513(33) 


0.9 


0.9(8) 


2.98(9) 


0.888(46) 


0.7(7) 


2.89(8) 


0.526(15) 


1.1 


0.8(7) 


2.97(9) 


0.881(41) 


0.9(7) 


2.85(8) 


0.509(21) 



ami = 0.006 





Dipole 


Tripole 


Q 2 cut [GeV^] 


XVdof 


F 2 (0) 


M~ 2 [GcV- 2 ] 


XVdof 


F 2 (0) 


My 2 [GeV~ 2 ] 


0.5 


1.7(1.5) 


3.14(7) 


0.797(39) 


1.6(1.5) 


3.10(7) 


0.492(22) 


0.6 


1.4(1.2) 


3.16(7) 


0.810(35) 


1.2(1.1) 


3.10(6) 


0.495(20) 


0.7 


1.5(1.1) 


3.18(7) 


0.825(33) 


1.1(1.0) 


3.12(6) 


0.501(19) 


0.9 


1.4(1.0) 


3.22(6) 


0.851(26) 


1.0(8) 


3.13(5) 


0.505(14) 


1.1 


1.3(9) 


3.24(5) 


0.861(22) 


1.0(7) 


3.11(5) 


0.499(12) 



ami = 0.008 





Dipole 


Tripole 


Q 2 cut [GeV^] 


XVdof 


F 2 (0) 


[GeV-] 


x7dof 


F a (0) 


M- 2 [GeV-' 2 ] 


0.5 


2.2(1.7) 


3.26(7) 


0.813(34) 


2.1(1.7) 


3.21(6) 


0.501(19) 


0.6 


1.6(1.3) 


3.26(6) 


0.813(33) 


1.7(1.3) 


3.20(6) 


0.497(19) 


0.7 


2.5(1.4) 


3.29(6) 


0.841(31) 


2.1(1.3) 


3.22(6) 


0.511(17) 


0.9 


2.1(1.2) 


3.31(5) 


0.851(22) 


1.8(1.1) 


3.22(5) 


0.506(12) 


1.1 


2.0(1.1) 


3.32(5) 


0.862(20) 


1.6(1.0) 


3.21(5) 


0.502(10) 
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Appendix B: SMEARED NUCLEON SOURCES FOR DOMAIN WALL FERMIONS 



Since careful optimization of the interpolating field for the nucleon source is crucial for the high precision calculations 
described in this work, in this appendix we describe in detail our optimization procedure and record the optimal 
parameters in two commonly used conventions. 

We have two objectives in constructing sources for propagators that will be optimal for calculating hadronic matrix 
elements. The first is to maximize the overlap between the interpolating field acting on the QCD vacuum and the 
hadronic ground state. The second is to minimize fluctuations arising from the source itself. Let N denote an 
interpolating field with the quantum numbers of the hadron, |^) = C _1 / 2 iV|ri) denote the normalized state obtained 
by its action on the vacuum, and \n) denote the n th eigenstate of the hadron (projected to zero momentum in the 
present discussion). Then, maximizing |(0| , I / )| 2 minimizes the contributions of excited states to the measurement of 
the hadronic matrix element of an operator O 

(N(t 3 )0(h)N(ti)) =C^(*|n)(n|0|?«)(m|*)e- £; ™( t3 - t2 )- is ™( t2 - tl ), (Bl) 

and hence enables one to reduce the source-sink separation while controlling contamination from excited states as 
discussed in Sect. IV Al 

The first objective is met by using smeared propagators and treating the rms radius of the smearing as a variational 
parameter. Although similar effects can be accomplished with gauge fixed sources, we use gauge invariant sources of 
the Wuppertal, or equivalently, Gaussian form by smearing a delta function source over the three spatial dimensions 
of the source time slice. 

Wuppertal smearing of a point source at the origin of time slice t is defined in the MIT USQCD software as 



i>(x, t) = 1 1 + a [U(x, i)S x+ly + U\x- i, i)S x _ ly 



N 



i = l 



(B2) 



and Gaussian smearing is defined in Chroma software as 



ip(x,t) 



o 2 V< 



N 



AN 

2S « 



ip(x,t) 



= 1- 



3(7 

2N 



a 2 /AN 
1 - 3cr 2 /2N 



N 



(B3) 



The Chroma and MIT parameters arc related by 



a 



a 2 /AN 
1 -ia 2 /2N 

2Na 
3a + 1/2 ' 



(B4) 
(B5) 



Note that there is an instability for a < 0, since the sign of the source generated in Eq. (|B2|) is then (—l) x+v+z , 
and the resulting spatially oscillating source has an extremely poor overlap with the physical ground state. In terms 
of the Gaussian parameters, the instability arises for N < 3cr 2 /2. 

Because the smeared sources contain link variables U, the statistical fluctuations in correlation functions using 
these sources are larger than those arising from point sources. To attain our second objective of minimizing the 
fluctuations arising from the source itself, it is highly advantageous to perform APE smearing of the gauge links used 
in generating the source on the time slice of the source. In each iteration of APE smearing, each link is replaced by 
a linear combination of itself and the sum of staples within that time slice, and projected back onto SU(3) as follows 



U 



(JV+l) 



= Proj 



SU(3) 



U 



N 



n tjn T T m 

x,j x+j,i x-\-i,j 



(B6) 
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and the APE smearing is repeated N times. An alternative notation is 



U 



(JV+l) 



Proj 



SU(3) 



ATT N + V/T^ TI N f/ wt 

^ x,i > / j ^ x,j^ x-\-j,i^ x-\-i,j 



so that 



A = 1//3. 

A convenient measure of the smearing of the source ip(x,t) in Eq. (|B2j) is the rms radius 



/ d 3 x\x\ 2 ip*(x, t)tp(x, t) 



J d 3 xip*(x, t)ip(x, t) 



1/2 



(B7) 



(B8) 



(B9) 



and Fig. 2 of Ref. [73[ shows how r rms depends on the parameters N and a. As one expects from the fact that smearing 
is a random walk governed by the gauge fields, the rms radius is approximately proportional to y/N. Since the size 
of the source is nearly independent of a for a > 3, at which point the constant term in Eq. (|B2|) becomes negligible 
relative to the hopping term, in all our calculations, we use a = 3, which provides the maximum r rms for a given 
number of smearing steps N. 

It is simplest to think about optimization criteria for Wilson fcrmions, for which one can construct a transfer matrix 
and correct propagators such that the two-point correlation function has quarks and antiquarks properly normal 
ordered at zero time separation [74]]. In this case, the source may be optimized straightforwardly by maximizing the 
overlap between the normalized state created by the action of the source \^^) = C~ l / 2 N^ where the source 

ATM 

has rms radius r, and the normalized ground state of the nucleoli |0). Denoting the momentum projected normalized 
eigenstates of the nucleon by \n) and their energies by E n , the momentum projected two-point correlation function 
may be expanded: 



C (r) (t) = J d 3 x{N {r \x,t)N {r \0,0))=C^\{¥^ 



■,-E-t 



(BIO) 



where C is an unknown normalization constant. Since one can directly measure the correlation function at zero time 
separation 



A (r) =C M( ) =CJ2\(* 



Ml 



(Bll) 



and reliably fit the large t behavior of the correlation function to extract the ground state contribution 

B (r) =C (frM|0) 2 , 
the probability that the source contains the nucleon ground state is given by 



5 (0 _ 



SM 



$M| ) 



(B12) 



(B13) 



Using this criterion, Bratt [75j| has recently shown that optimizing the source size with 4-component nucleon sources 
yields a maximum overlap of 35%, projecting onto the upper two components (in the Bjorken-Drell convention for 
which these components yield the non-relativistic limit) increases the overlap to 50%, and APE smearing of the gauge 
links in the source further increases the overlap to 80%. 

For domain wall fcrmions, which do not have a local transfer matrix, we consider the following generalization of 
Eqs. (|B1HIB13|) . which compares the ratio of the correlation function and extrapolated ground state contribution at 
time t instead of time 0: 



A {r) (t) 
B^ r) {t) 

r {r) (t) 



C (r \t), 
C (* M |0) 



-E t 



(B14) 
(B15) 

(B16) 
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Figure 22: The rms radius of a gauge invariant smeared source as a function of the coefficient a and number of smearing steps 
N defined in Eq. ()B2|I for the coarse (left panel) and fine (right panel) lattices. The curves projected in the horizontal plane 
show the numbers of smearing steps required for the specific values of rms radii shown in the key. 



This ratio, V^ r '(t), ranges from the overlap at t = to 1 in the limit t — > oo. We expect that for small t, it is 
still a good measure of the presence of excited state components in the source and should have a maximum close to 
the maximum in V^ r \ This expectation is borne out in the case of Wilson fcrmions, and we note that this criterion 
gets even better as the lattice spacing decreases. Since we are only interested in the dependence of V^ r '(t) on the rms 
radius r and the absolute normalization for i ^ has no physical significance, it suffices to calculate the following 
ratio for large to 

For each value of t, it is convenient to normalize the curve such that its maximum value is unity. Hence, defining the 
rms radius at the maximun as r* , our final criterion for optimizing the smearing is the ratio 

R {t) ~ C(r*)(to)/C(r-)(ty (B18) 

Equation (|B18|) has the computational advantages that all oscillating terms in the time dependence of the correlation 
functions cancel out of the ratios and that jacknife or bootstrap analysis enables accurate measurements on small 
ensembles. 

We now show the results of optimizing the ratio R^ r '(t) on a coarse 24 3 x 64 domain wall lattice with a — 0.114 fm, 
TOjr = 420 MeV, m s — 0.04, and m u — 0.01, using 32 configurations and on a fine 32 3 x 64 domain wall lattice with 
a = 0.081 fm, = 310 MeV, m s = 0.03, and m u = 0.004 using 33 configurations. We included both APE smearing, 
with (3 = 0.3509, and Wuppertal smearing with a = 3. Because APE smearing smoothes the links, the rms radius 
obtained from a given number of Wuppertal steps changes with the number of APE steps, becoming slightly larger 
as the number of APE smears increases. Figure [55] shows the rms radius calculated as a function of both the number 
of APE and Wuppertal steps for both lattice spacings. 

Since fluctuations in the normalization of the source directly contribute to the overall fluctuations in correlation 
functions, it is desirable to use APE smearing to smooth the spatial links used in generating the source and thereby 

diminish the fluctuations. A simple measure of these fluctuations is the relative fluctuation = /p? , where 
O is the rms radius defined in Eq. (|B9[) . Figure 1531 shows the dramatic effect that APE smearing has in reducing these 
fluctuations for both lattice spacings. Since the incremental benefit of successive smearing becomes small beyond 25 
smearing steps, we have chosen to use 25 steps throughout. Note that for the largest number of Wuppertal steps, this 
reduces the noise by a factor of more than 5 in each case. 

Figure [53] shows the primary result of the calculation for both lattice spacings. For the coarse lattice, the ratio 
R( r >(t) is calculated at six values of the number of Wuppertal steps, N = 10,20,30,50,70,100, corresponding to 
rms radii, r = 2.07, 2.89, 3.51, 4.46, 5.19, and 6.06 lattice units respectively. We chose r* = 4.46 fm and calculated 
bootstrap error bars using 32 configurations. Instead of normalizing at a single value of to as in Eq. (p!8jl . the errors 
in the ratios in Fig. [5H were further reduced by normalizing to an exponential fit to each correlation function in the 
region t = [6—12]. These results are completely consistent with those of a single to, but display the shape of the 
maxima more precisely. Note that for all four values £ = 1,2, 3, and 4, the curves are accurately determined and the 
ratio R( r \t) has a maximum at approximately the same point, r = 4.0, corresponding to N — 40. Thus, we believe 
our optimization criterion is robust and statistically accurate for domain wall fcrmions. 
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Figure 23: The source variance SO/O as a function of the rms radius and number of APE smears N for the coarse (left panel) 
and fine (right panel) lattices. The curves projected in the horizontal plane show the numbers of smearing steps required for 
the specific values of source variance shown in the key. 




2 2 

sqrt(<r >) sqrt(<r >) 



Figure 24: The left panel shows the ratio Iv r \t) for t = 1,2,3, and 4 as a function of r on the coarse lattice. The solid curves 
are splines passing through the mean values to guide the eye. This graph provides a robust determination of the optimal rms 
radius r = 4.0 lattice units, corresponding to N = 40 Wuppertal smearing steps. The right panel shows the analogous ratio 
R( r \t) for t = 1 and 2 as a function of r on the fine lattice. 

For the fine lattices, the ratio R( r \t) is calculated at 5 values of the number of Wuppertal steps, N = 30, 50, 70, 100, 
and 150, corresponding to rms radii r = 3.76, 4.77, 5.56, 6.51, and 7.77 lattice units respectively. We chose r* = 5.56, 
normalized by exponential fits to each correlation function in the region t — [6 : 12], calculated jackknife error bars, 
and only included t = 1 and 2 to avoid making the graph confusing due to the larger error bars. The maximum 
occurs at approximately r = 6.0 lattice units, corresponding to 84 Wuppertal smearing steps. This result appears 
reasonable, since assuming a constant rms radius in physical units would imply that the the rms radius on the coarser 
lattice of 4.0 lattice units would scale to 4 x 0.123/0.093 = 5.3 lattice units on the present lattice, and the pion mass 
on the finer lattice is somewhat lighter. 

We summarize the final parameters for optimal sources used in this work in Tab. IXVIII1 where the parameters are 
defined in Eqs. [[B31- |B8]). 



lattice 


APE 


smearing 


Wuppertal smearing 


size 


a (fm) 


P 


A 


Nape 


Q 


a 


N w 


(r ) 


0.114 


0.3509 


2.85 


25 


3 


5.026 


40 


4.0 


0.084 


0.3509 


2.85 


25 


3 


7.284 


84 


6.0 



Table XVIII: Parameters for optimal sources. 
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